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Abstract In this chapter, we describe three related studies of the universal physics of 
two-component unitary Fermi gases with resonant short-ranged interactions. First we 
discuss an ab initio auxiliary field quantum Monte Carlo technique for calculating 
thermodynamic properties of the unitary gas from first principles. We then describe 
in detail a Density Functional Theory (DFT) fit to these thermodynamic properties: 
the Superfiuid Local Density Approximation (SLDA) and its Asymmetric (ASLDA) 
generalization. We present several applications, including vortex structure, trapped 
systems, and a supersolid Larkin-Ovchinnikov (fflo/loff) state. Finally, we discuss 
the time-dependent extension to the density functional (TDDFT) which can describe 
quantum dynamics in these systems, including non-adiabatic evolution, superfiuid to 
normal transitions and other modes not accessible in traditional frameworks such as 
a Landau-Ginzburg, Gross-Pitaevskii, or quantum hydrodynamics. 
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9.1 Introduction 

The question of how pairing correlations between two types of fermions develop 
with interaction strength has fascinated physicists for decades, beginning with the 
papers of Eagles [ 1 1 and Leggett |2 1, and followed by many others |[3][4][5][6]]. These 
initial studies focused on the inter-species pairing gap at various temperatures as the 
pairing interaction varied throughout the entire BCS-BEC crossover from weak to 
strong attraction. 

Eagles and Leggett AUG]! solved the Bardeen-Cooper-Schrieffer (BCS) mean-field 
equations only in the particle-particle (pairing) channel: The prevailing attitude 
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(influenced by electronic systems) was that the pairing gap is much smaller than the 
self-energy (exponentially suppressed in weak-coupling), and that the presence or 
absence of pairing correlations was a tiny effect compared to the background density 
which determined the self-energy. Subsequent improvements to the theory focused 
only on a more accurate description of the pairing channel Il3ll4l I5ll6l 171 [8l [9l ITOl ITTI . 
neglecting the so called "Hartree-Fock" contributions to the total energy of such a 
system. 

However, even in the weak coupling limit (a < and kp\a\ <C 1 where the Fermi 
energy £p — p 2 F /2m, the Fermi momentum pp = hkp =h(37t 2 n) 1 ^, n is the total den- 
sity, and a is the two-body s-wave scattering length) — which was rather thoroughly 
studied in the 1950's |[T2l[T3l — it was evident that the "Hartree-Fock" and higher 
order particle-hole contributions dominate in the total energy. These contributions 
can be described perturbatively in terms of the small parameter kpa (in both BCS and 
BEC limits). In the BCS limit, for example, the leading contributions enter at linear 
order °< Epkpa while the particle-particle pairing contributions are exponentially 
suppressed °= Ef exp(n/kpa). 

Despite neglecting the dominant particle-hole contributions, these mean-field 
studies correctly captured many of the qualitative features of the BCS-BEC crossover. 
This can be partially attributed to the fact that the particle-particle channel correctly 
accounts for the two-body bound state that dominates in the extreme BEC limit 
at strong attraction (however, higher order effects — describing the dimer-dimer 
interaction for example — are not correct lfT^[T5l[l6l[T7l ). 

At unitarity, the majority of the interaction energy is due to the particle-hole 
channel: see [TTsTl where the energy of the normal state at T = was evaluated 



for the first time and the discussion in Sec. 9.2.6 In particular — above the critical 
temperature T c , for example — the total energy of the normal phase exceeds the ground 
state energy by only about 20% or so fl9l : This means that the condensation energy 
gained by the particle-particle pairing interaction is a relatively small contribution to 
the total interaction energy. A quantitative description of unitary physics must thus 
include these "Hartree-Fock" contributions and go beyond the simple mean-field 
models used initially to study the crossover. 

In 1999, G. F. Bertsch [20| emphasized the special role played by the problem 
of a two-species Fermi gas at unitarity with large scattering length. In the crust of 
neutron stars one can find a situation where the scattering length a of the interaction 
is anomalously large compared to the other length scales, the average interparticle 
separation n~ 1//3 , and the range ro of the interaction: ro <C n~ 1,/3 <C \a\. Since the 
Fermi momentum is small (kp ro <C 1), the neutrons effectively interact only in the 
relative s-partial wave, and the ground state energy should be some function of 
the physical parameters defining the system E gs — f(N,V,ro,a,h,m), where N is 
the particle number contained in a volume V of the system. In the formal limit of 
kFfo — > and 1 jkpa — > this function simplifies: 

E gs =f(N,V,h,m) = ^£ F N$, (9.1) 



3 



9.1 - Introduction 



and all the non-perturbative effects are described by a single dimensionless constant: 
4 (often referred to as the Bertsch parameter). At finite temperatures the total energy 
of the system becomes a slightly more complicated function, since now it depends 
also on the temperature T: 

E(T)=f(N,V,T,h,m) = 3 -e F N^ (J-^j . (9.2) 

The Bertsch parameter (along with all other thermodynamic properties) becomes a 
"universal" function of the dimensionless variable T/ef ET|. 

In 1999 it was not yet clear whether this limit existed: One might expect such 
a system to collapse, since the naive coupling constant g = Anh 2 a/m is infinite at 
unitarity. Baker 11221 l23l l24l provided the first clue that this system was actually 
stable. Carlson and collaborators [18| subsequently calculated the energy of this 
system, proving that it was stable, and showing that the superfluid paring gap was 
very large. Meanwhile, using a Feshbach resonance to induce an extremely large 
scattering length, J. E. Thomas and his collaborators 11251 produce for the first time a 
quantum degenerate unitary gas of cold-atoms in a trap, thus providing experimental 
evidence that this system is indeed stable. 

There has since been an explosion in both theoretical and experimental studies of 
resonant Fermi gases near the unitary regime (see for example the reviews [26 27 28 
|29l[30)). On one hand, cold-atom experiments can simulate other systems of interest; 
for example, dilute superfluid neutron matter which can only exist in the crust of 
neutron stars, various condensed matter systems (the unitary gas exhibits a pseudogap 
that might shed light on the pseudogap in high-temperature superconductors), and 
quantum systems with extremely low viscosity similar to quark-gluon plasmas 
observed in ultra-relativistic heavy-ion collisions. On the other hand, the simplicity 
of the system provides an excellent vehicle through which the plethora of many-body 
techniques can be put to rigorous test, including both traditional approaches, as well as 
modern developments such as the e-expansion and Ads/CFT correspondence ||3Tl[32l , 

We shall not provide a cursory review of current theoretical techniques, but will 
instead focus on a couple of theoretical methods that have produced a large reliable 
set of information about the properties of unitary Fermi gases. The first approach 
is an ab initio Quantum Monte Carlo (QMC) method that has accurately evaluated 
many properties of these systems, and has been confirmed by experiments. The 
second approach is Density Functional Theory (DFT), which is in principle, an exact 
approach commonly used for describing "normal" systems (no superfluidity). We 
show how to extend the DFT to describe both superfluid systems and time-dependent 
phenomena, and how the DFT allows us to address phenomena that are essentially 
impossible to describe within a QMC approach. 
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9.2 The Quantum Monte Carlo Approach 

9.2.1 From the Physical Problem to the Lattice Formulation 

Atomic collisions in a trap occur at very low relative velocities (due to the diluteness 
of the gas) and this fact allows us to restrict the description to using the lowest 
partial waves only. In practice, the s-wave scattering phase shift fully determines the 
properties of a unitary Fermi gas, for which ro -C «~'/ 3 <C \a\. The detailed physics 
of the collision is more complicated since atoms are not point-like objects and can 
appear in various configurations. Roughly speaking, these can be associated with 
various valence electronic configurations. For example, atoms with a single valence 
electron (such as 6 Li) can form two possible electronic configurations in a binary 
system: a singlet and a triplet configuration. The inter-atomic potential describing a 
singlet configuration corresponds to the symmetric spatial wave function. It admits 
the existence of a bound state and corresponds to the closed (inaccessible) scattering 
channel. The triplet channel, on the other hand, is open and shallow: due to large 
(mainly electronic) magnetic moment, its energy can be easily tuned with respect 
to the closed singlet channel threshold by adjusting an external magnetic field. This 
allows experimentalists to use a Feshbach resonance to tune the effective interaction 
in the open channel to virtually any value: in particular, experiments with dilute 
clouds of cold atoms can directly probe the unitary regime. 

A typical Hamiltonian describing the two channel atom-atom collision is of the 
form 051 134ll35ll36ll37l: 

H =^ + t i V l f + V >) + VC + V "> < 9 - 3 ) 

where M r is the reduced mass of two atoms, V h f — a^/S^ ■ S n /h 2 is a hyperfine 
interaction term for each atom (with hyperfine constant a/,/0, and S e and S" are 
the total electron spin and the total nuclear spin respectively. The Zeeman term 
V z = (YeSf + Y n S'l)B describes the interaction with the external magnetic field B 
which is assumed to be parallel to the z-axis. The terms V c and V d denote the 
Coulomb interaction and dipole-dipole magnetic interaction, respectively. The dipole 
term contributes weakly to the interaction and can be neglected. The Coulomb term 
distinguishes singlet and triplet channels (due to different symmetry properties of 
electronic wave function) and produces different interaction potentials in both chan- 
nels. Consequently, the continuum of the singlet channel lies above the continuum of 
the incident triplet channel. At very low collision energies, only the singlet channel 
is open. However the hyperfine interaction couples the singlet and triplet states and 
consequently, resonant scattering may occur due to the bound state of the singlet 
potential (see reviews Il38ll39l l40l and references therein). An external magnetic field 
can thus be used as an experimental knob to control the resonance position, effec- 
tively altering the atom-atom collision cross-section. In the limit of low collisional 
energy, the effective scattering length for two colliding atoms is well described by 
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a(B)=a + C , (9.4) 

a — Ores 

where «o is the triplet channel off-resonant background scattering length, and C > 0. 
The second term results from the coupling to the closed channel, and B les is the value 
of the magnetic field where the Feshbach resonance occurs. In this way, experiments 
may realized the unitary Fer mi g as by considering dilute systems (ro <C n~ 1,/3 ) and 



tuning the scattering length (9.4 1 near the resonance (n -1 / 3 <C \a\). 

To determine the thermodynamic properties of an ensemble of fermionic atoms in 
a non-perturbative manner, we consider the system on a three dimensional (3D) cubic 
spatial lattice with periodic boundary conditions. The system consists of two species 
of fermions that we shall denote "a" and "b". In dilute neutron matter these would 
correspond to the two spin states of the neutrons, while in cold atom experiments 
these are the two populated hyperfine states. Although there are physical processes 
that can convert one species to another, for the purposes of the experiments we shall 
describe, these transitions are highly suppressed and one can consider each species 
to be independently conserved. 

The lattice spacing I and size L — N s l introduce natural ultraviolet (uv) and 
infrared (IR) momentum cut-offs given by hk c = nh /I and HAq = 2nh /L, respectively. 
The momentum space has the shape of a cubic lattice, but in order to simplify the 
analysis, we place a spherically symmetric UV cut-off, including only momenta 
satisfying k < k c < n/l. In order to minimize the discretization errors, the absolute 
value of scattering length must be much larger than the lattice spacing: a^> I. 



9.2.2 Effective Hamiltonian 

As discussed in the introduction, it has by now been well established that the unitary 
regime exists and is stable. Hence, any sufficiently short-ranged interaction with 
large scattering length will exhibit the same universal physics. Here we use a contact 
(zero-range) interaction V{r\ — T2) = — g8(ri — T2) regularized by the lattice, which 
defines a momentum cut-off hk c . (We require all two-body matrix elements to vanish 
if the relative momentum of the incoming particles exceeds this cutoff.) The second 
quantized Hamiltonian of this system is 

H = Jd\(^- £^ f+(r)^%(r) +g« a (r)%(r)^ , (9.5) 

where n a (r) = i//+(r)y/<j(r). Once the cutoff is imposed, the value of the bare 
coupling g can be tuned to fix the value of the renormalized physical coupling — in 
this case, the s-wave scattering length a. The relation between a and the coupling 
constant g can be obtained from T matrix describing two-particle scattering induced 



by the interaction (9.5 1 with the s-wave phase shift: 
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Anh 2 2 k 

kcoto = k c In 

gm n n 



k c ~k 



k c + k 



The low-momentum expansion of the scattering amplitude reads: 



(9.6) 



/(*)< 



-ik- 



Antf 
gm 



2k c 2k 2 

K Kk c 



0(k 3 ) 



(9.7) 



At low momentum we have f(k) = [— ik— l/a + r e ffk 2 /2 + 0(£ 3 )]~ 1 , which gives 
the relation between the bare coupling constant g and the scattering length a at a 
given momentum cutoff hk c : 



1 m k c m 



g Antra 2n 2 tr Antra 



2k c a 



(9.8) 



One has to remember, however, that the value of the coupling constant g has been 
determined for the two body system in its center of mass frame. On the other hand 



the Hamiltonian (9.5 i is supposed to describe an ensemble of fermions in the box. 
Consequently, only a fraction of interacting pairs have their center of mass at rest 
with respect to the box. Most of the interaction processes will occur for pairs for 
which the center of mass velocity is nonzero. It implies that their mutual interaction 
will be characterized by a slightly different scattering length than l |9.8| >. Consequently, 
the Hamiltonian will generate a systematic error in the description of interacting 
fermions. This error will scale as kf/k c and in order to minimize its influence one 
should keep the particle density as small as possible. Another source of systematic 
error is related to the nonzero effective range, which is generated by the interaction 
and is independent of the coupling constant r e ff = A/(nk c ). Note however that the 
choice of k c described above implies that r e g- < I. 



9.2.3 The Hubbard-Stratonovich Transformation 

Since we are interested in the finite temperature thermodynamic properties of the 
system, it is natural to use the grand canonical ensemble to evaluate physical quanti- 
ties. This is equivalent to considering a small portion of volume V = L? in thermal 
and chemical equilibrium with the larger system. Consequently we allow for energy 
and particle exchange between our subsystem and the larger system, fixing only 
the average values of these quantities in the box. The thermodynamic variables are 
thus the temperature T, the chemical potential ji, and the volume V. The partition 
function and average of an observable O are calculated according to 

Z(p,ii,V)=TT{exp[-P(6-iift)]}, 
s Tr {<3exp[-B(#-U#)]j 
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where J3 = 1/T (in this work we will take Boltzmann's constant to be kg = 1 so that 
temperature is expressed in units of energy). In order to be able to calculate these 
quantities we first factorize the statistical weight using the Trotter formula: 



exp[-P(H-nN)] = IJexp[-T(i7-jU#)] 



(9.10) 



where j3 = A^T. The next step is to decompose the exponentials on the right hand side 
into exponentials that depend separately on the kinetic and potential energy operators. 
The second order expansion is (higher orders require more effort, see flTll42ll43ll44l '): 



exp[— — fiff)] 



- exp 



~2 



exp(— TV)exp 



2 



0(t 3 ), (9.11) 



where K is the kinetic energy operator, whose dispersion relation, for momenta 
smaller than the cut-off, is given by = h 2 k 2 /2m. Since T has the dimension of in- 
verse energy, the above approximate representation makes sense only if T max \\V\\ <C 1 
and T max 1 1 K — flN\\ <C 1. Since both the interaction and kinetic energies are extensive 
quantities, this restriction might appear as very strict. However, after performing a 
Hubbard-Stratonovich transformation (see below), this restriction is considerably 
eased and both the kinetic and the interaction energies in these inequalities are re- 
placed by the corresponding intensive energies per particle. It is important to note 
that, because we have used the expansion up to 0(t 3 ), when calculating the partition 
function the error becomes 0(t 2 ). Indeed, the statistical weight involves a product 
of N T factors and is given by the following expression: 



exp[-/3 (H - nN)} = exp 



x(K-nN) 



(Si . \ 

x |Jexp[— T^]exp[— t(£ — pf})] exp 



t(K-hN) 



-0(t 2 ) (9.12) 



Note also that this approach does not depend on the choice of dispersion relation in 
the kinetic energy term. However various choices of representation of derivatives 
on the lattice may lead to different discretization errors B31 . In our case we shall 
consider the kinetic energy operator in momentum space, e(k) — h 2 k 2 /2m, which 
minimizes the discretization errors. 

In order to efficiently evaluate the term containing the interaction, one has to 
replace it by the sum (or integral) of one body terms. This can be done with the 
Hubbard-Stratonovich transformation [46 1 . The transformation is not unique, and 
we take advantage of this freedom to ensure an efficient summation (or integration) 
scheme. In our case, due to the simplicity of the interaction term, a discrete Hubbard- 
Stratonovich transformation can be applied, similar to that in |47|: 
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exp[-gT« fl (r)n 6 (r)] = - £ [1 +Aa(r,T/)n a (r)][l +Aa(r,T/)^(r)], 



(9.13) 

where A = ^/exp(— gr) — 1, T,- labels the location on the imaginary time axis, j = 
1 , . . . ,7V T , and ci(r, T;) is a field that can take values ±1 at each point on the space- 
time lattice. This identity can be proved simply by evaluating both sides at n^ a b j (r) = 
0, 1. This discrete Hubbard-Stratonovich transformation is sensible only for A < 1, 
which means that the imaginary time step cannot exceed log2. The advantages 
of this transform is discussed, for example, in Il47ll45l . 

Taking all this into account, the grand canonical partition function becomes 



Z(p,n,V) = Tr{eKp[-p@-nf))]}= Yl®<j(r,Tj)TvW({a}), (9.14) 
where we define 



A'r 



^(M)=n^(M) 

7=1 



and 



^•({a})=exp 





t(K-iiN)' 


J exp 


2 



Since a(r, t) is discrete, the integration is in fact a summation: 



/n^^i— e e 



)}=±1 {<7(r,T,v r )}=±l 



where 



E E E E 

{a(r,T/)}=±l cj((1,0,0),t / )=±1ct((2,0,0),t / )=±1 o((N s ,N s ,N s ),Tj)=±l 



dz[h({(j})-nN]}, 



(9.15) 



(9.16) 



(9.17) 



(9.18) 



In a shorthand notation we will write 
<#({cr})=T T exp 



where T T stands for an imaginary time ordering operator and h({a}) is a resulting a- 
dependent one-body Hamiltonian. It is crucial to note that ^ ({d}) can be expressed 
as a product of two operators which describe the imaginary time evolution of two 
species of fermions: 
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W({o}) = %({o})W a ({o}), (9.19a) 
^(M) = ri^(M), €(M) = n# ;£! ({(T}). (9.19b) 

As we only consider unpolarized systems, for which jj, a = jib — jj, , the operators for 
both species a and b are identical. 

The expectation values of operators take the form: 

> z{p,n,v) 



nij@e(r,Tj)TrW({cj}) Tr(M({cr}) 
Z(fi,il,V) Tr^({cr}) : 



(9.20) 



where we have introduced Tr fy({a}) for convenience: in the numerator it represents 
the probability measure used in our simulations (see below), and in the denominator 
it serves the purpose of moderating the variations of Tr O'ftf ({cr}) as a function of 
the auxiliary field a. 

All of the above traces over Fock space acquire very simple forms [48 49), and 
can be easily evaluated. In particular, Tr ^ ({c}) can be written as 

Tr#({cr})=det[l + ^({c7})] =det[l + ^({a})]det[l + i %({cr})], (9.21) 

where ^/ (without the hat) is the representation of tyt in the single-particle Hilbert 



space. The second equality is a result of the decomposition (9.19 1 and is easy to 
prove by expanding both sides. For symmetric (unpolarized) systems the chemical 
potentials \i a = are the same for both species of fermion, so it follows that 
det[l +%({a})} = det[l +<%({cr})]. This implies thatTr #({cr}) is positive, i.e., 
that there is no fermion sign problem. Indeed, this allows to define a positive definite 
probability measure: 

_ Tr#({cr}) _ {det[l + ^,({(7})]} 2 



z(p,n,v) z(p,n,v) 
1 



exp(2tr(log[l +%({cj})])) (9.22) 



where the exponent in the last equation defines the negative of the so-called effective 
action. The positive definite probability measure is crucial for Monte Carlo (MC) 
treatment, allowing for statistical sampling of the a space. When considering the 
polarized system, the sign problem inevitably occurs, making the Monte Carlo 
procedure very difficult. The sign problem appears also when more complicated 
forms of interaction are applied. In such a case one can sometimes cure the problem 
by properly choosing the Hubbard-Stratonovich transformation |50|. 

The many-fermion problem is thus reduced to an Auxiliary Field Quantum Monte 
Carlo problem (AFQMC), to which the standard Metropolis algorithm can be applied, 



10 



9.2.3 - The Hubbard-Stratonovich Transformation 



using ( |9.22[ ) as a probability measure. Before moving on to the details of our Monte 
Carlo algorithm, we briefly discuss the expressions used to compute a few specific 
thermal averages. 

Let us consider the one body operator 



0= E /dVitf 3 r 2 ft+(ri)0„(ri,r 2 )w(r 2 ) 



i.t=b,c 



(9.23) 



From ( |9.20| i it follows that 

{0) = £ } P({ * }) ^({^ = { ? } P({ff}) det[l + ^(M)]- 
The calculation of the last term requires the evaluation of 

Tr[^+(n)^(r 2 )#({a})] - 5 4 ,det[l + ^({a})] 2 « 4 (r,,r 2 ,{(j}) 
where s and t run over both species (a or b), and 



n s (n,r 2 ,{a}) = E <iM r i) 



i+^(M) 



(9.24) 



(9.25) 



(9.26) 



k,,k, 



Here <Pk(r) = exp(/k • r)/L 3 / 2 are the single-particle orbitals on the lattice with 
periodic boundary conditions, and hence quantized momenta k = 2nn/L. This holds 
for any 1-body operator O, if % is a product of exponentials of 1-body operators, 
as is the case once the Hubbard-Stratonovich transformation is performed. It is then 
obvious that the momentum representation of the one-body density matrix has the 
form 

%{{<*}) 



« s (ki,k 2 ,{cr}) = 



\ + %{{a}) 



(9.27) 



ki.ko 



which, for a non-interacting homogeneous Fermi gas, is diagonal and equal to the 
occupation number probability 1 /(exp[j3(£k — fl)] + 1) of a state with the energy 
E k = h 2 k 2 /(2m). 

Summarizing, the expectation value of any one-body operator may be calculated 
by summing over samples of the auxiliary field c(r, T/): 



(O) 



Y[9G{Y,Xj)P{{a}) £ £ 0„(ri,ri)f!,(ri,r 2> {a}) (9.28) 



In particular, the kinetic energy can be calculated according to: 



11 



9.2 - The Quantum Monte Carlo Approach 



n r , Tj ^cj(r,T ; )Tr^({a}) TrKW{{o}) 



Z(P,n,V) 



Tr^({cr}) 

k<k c 



= /n^(r,T,)P({cr})£ r £ 

•J r T • t « 



n s (k,k,{(7}) 



2m 



(9.29) 



(9.30) 



Analogously, for a generic two-body operator: 
6 = E / dV;dV' 2 dV 1 dV 2 ^+(r;)^+(r , 2 )O s( „ v (r;,r , 2 ,r 1 ,r 2 )^(r 2 )v/„(r 1 ; 

In order to calculate (O) one needs to evaluate the expression 
Tr^+Cri)^^)^^)^^)^^^})] 

= (det[l + ^({cr})]) 2 ^5 SH ^ v n,(ri,n,{a}H(r2,r2,{a}) 

- 8 sv 8 tu n s (T f 1 ,r 2 ,{o})n t (T' 2 ,r 1 ,{o}) S J . (9.31) 
Hence, for the expectation value of the two body operator we get 
(0)= [H®o(r,Tj)P({o}) 



x E E 

rpr^ri ,r 2 s,t=a,b 



0^(r' 1 ,r 2 ,ri,r 2 )n s (r' 1 ,ri,{(7})nf(r 2 ,r2,{(7}) 



- O sm (r\ , r' 2 , n , r 2 )n. s (r' 1 , r 2 , {a})n t (r 2 , n , {a}) 
In particular, the expectation value of the interaction energy reads: 

= -* / n^ a ( r '' r 7) p ({^})E' i «( r ' r '{ (j }H( r ' r '{ <7 }) 

It should be noted that in the symmetric system (p, a = jj,b) 
n«(r,r',{a}) =n b {r,r' ,{a}). 

Hence, 

(V) = ~8 [ n^(r 1 T i )P({a})£K(r 1 r>})f 

J r,Xj r 

It is useful to introduce the correlation function 



. (9.32) 

(9.33) 

(9.34) 
(9.35) 
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8i(r) = y-J J d 3 ri d 3 r 2 (v/i(ri + r)v/ 6 t (r 2 +r)%(r 2 )i// a (ri)) 

=(|)7n^ ci(r '^ )jP({cj}) 

x j A 3 r l A\ 2 na{ri+r 1 r ll {a})n h {r 2 +r.r2.,{a}) 1 (9.36) 

(where N is the average particle number) which is normalized in such a way that for 
a non-interacting homogeneous Fermi gas g 2 (r) = 3j\ (kpr) / (kpr) and g 2 (0) = 1 . 



9.2.4 Stabilization of the Algorithm for Small Temperatures 



Once we have written the observables as in (9.20 1, the next step is to sum over all 
possible configurations of c(r, T/). This is still an impossible task, as for example, a 
lattice size N% x N T (where typically N x = 8 and N t ~ 1000), requires performing 
the sum over the 2 * xNt points in configuration space. It is in these cases that a 
Monte Carlo approach becomes essential. By generating JV independent samples of 



the field c(r, T/) with probability given by (9.22 1, and adding up the values of the 
integrand at those samples, one can estimate averages of observables with 0(1/V JV) 
accuracy. 

The standard Metropolis algorithm is used to generate the samples. Namely, at 
every MC step, the sign of a is changed at random locations of the space-time lattice 
(see lfT9ll45ll5TI for details). This procedure allows to probe the sigma space, in order 
to collect the set of statistically uncorrected samples. 

In order to compute the probability of a given a configuration, it is necessary 
to find the matrix elements of °i/ ', which entails applying it to a complete set of 
single-particle wave-functions. For the latter we chose plane waves (with momenta 
hk < hk c ). This choice is particularly convenient because one can compute the overlap 
of any given function with the whole basis of plane waves by performing a single 
Fast Fourier Transform (FFT) on that function [4-5 1 . 

The procedure described above requires many matrix multiplications to calculate 
°i/ . In particular at low temperatures the number of matrix multiplications grows 
rapidly and the matrices have elements that vary over a large range of magnitudes. To 
avoid numerical instabilities it is necessary to separate the scales when multiplying 
the matrices, and a more costly but robust algorithm such as the Singular Value 
Decomposition (SVD) is required. In this section we follow the same approach 
developed in fl48l to introduce the SVD to our calculations. 

Let us write the matrix a i/ ({a }) more explicitly: 

^(M) = fl ^/(M) = ' ' ' «i , (9.37) 

7=1 
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where the #/t({<7}) are N x N matrices, for a single-particle basis of dimension N. 
Let us then define 

<% = 1 (9.38) 



To separate the scales one decomposes the matrix °i/ n -\ before multiplying it by 
W n to get This process begins as follows 

<&b = 1 (9.39) 

where and V\ are orthogonal matrices (not necessarily inverses of each other), and 
@i is a diagonal positive matrix containing the singular values of . The idea is that 
the actual multiplication should be done by first computing the factor in parenthesis 
in the last equation. This factor is then decomposed into ySli^i, in preparation for 
the multiplication by #3, and so on. A generic step in this process looks like: 

%> = Wn-1 - W n .?n-\%-\ Vn-\%-2 (9-40) 

so that in the end 

w Nx = w{{o}) = y N ^r Nx y Nx - x — v x = &&y, (9.41) 

where we have decomposed the full product in the last step. Calculating the determi- 
nant, and therefore of the probability measure, is straightforward if we perform one 
final more SVD in the following chain of identities: 

det(l + <2f ({a})) = det(l + 2>®V) = det(^(^V t + ®)V) 

= det(y yg>f y) = det(y y) det(i>) det(r f) (9.42) 

For equal densities (the symmetric case) we need this determinant squared, so we 
only care about the factor in the middle of the last expression: the other two factors 
have unit magnitude. Indeed, in that case we can write the probability measure as 

P({cr})=exp^f log^ (9.43) 

where dj > are the elements in the diagonal of §1, and M is the dimension of 
the single particle Hilbert space. The number of SVD's required to stabilize the 
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calculation grows as we increase j3. In our calculations we have made limited use of 
the SVD, ranging from 2 decompositions at the highest T to 8 decompositions at low 
T's. 



9.2.5 Finite Size Scaling 

The Monte Carlo calculations are performed in a box of finite size with a finite 
average number of particles. We are interested, however, in the thermodynamic limit 
N — > °°,V — > °° and N/V = const, so we need to consider the finite size scaling 
of the system so we can properly relate the values calculated in the box to their 
thermodynamic counterparts. This becomes particularly important in the vicinity 
of phase transitions where the correlation length <i; corr characterizing the non-local 
degree of correlation of a system diverges: 

^T-kr, (9.44) 

where t = 1 — T/T c , T c is the critical temperature, and v is a universal critical expo- 
nent. For the U(l) universality class, (which contains superfluid phase transitions), 
this exponent is well-known: v = 0.671. 

When dealing with systems that have a finite size L 3 , the theory of the renormal- 
ization group (RG) predicts a very specific behavior for the correlation functions 
close enough to the transition temperature (see e.g. 11521 ). In particular, the two-body 
density matrix K(L, T) that gives the order parameter for off-diagonal long-range 
order, scales as 

R(L, T) = L 1+ ^K(L, T) = f(x)(l+cL- m + •••), (9.45) 

where rj = 0.038 is another universal critical exponent, f(x) is a universal analytic 
function, x = (L/ t, corr ) l l v , and c is a non-universal constant, and co ~ 0.8 is the 
critical exponent of the leading irrelevant field. One should keep in mind that typically 
one knows neither c nor T c , but is interested in finding the latter. 

In a typical Monte Carlo calculation K(L, T) is computed for various lengths 
Li and temperatures T. The procedure to locate the critical point (characterized 
by scale invariance) involves finding the "crossing" temperatures Tn, for which 
R(Lj, Tij) — R(Lj,Tij) at two given lengths L, and Lj. Assuming that one is close to 
the transition (so that the correlation length is large compared to any other scale), 
one can expand /(jc(|*|)) = /(0) + f'{0)L x l v b\t\ (where we set t, corr = b\t\~ v ), and 
derive the relation 

iTc-Ti^KgiL^Lj), (9.46) 



where 

g (L i>L j)=L-:^ 



1 - 



l/v 



(9.47) 
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and K = cT c f(Q)/bf'(0). If there were no non-universal corrections to scaling (i.e. if 
c = 0), then fc = and T c = Tfj, which means that, upon scaling by the appropriate 
factor (as above) all the curves K(L, T) corresponding to different L's would cross 
exactly at T c . In general these corrections are present, and it is therefore necessary 
to perform a linear fit of 7y vs. g(Li,Lj) and extrapolate to infinite L in order to 
determine the true T c [45 1. 



9.2.6 Results: the Energy and the Entropy 

The results of our Monte Carlo simulations are shown in Figs. |9.l| and |9.2| [ 19, 
|5Tll45l . The Monte Carlo autocorrelation length was estimated (by computing the 
autocorrelation function of the total energy) to be approximately 200 Metropolis 
steps at T ~ 0.2ep . Therefore, the statistical errors are of the order of the size of the 
symbols in the figure. The chemical potential was chosen so as to have a total of 
about 45 particles for the 8 3 lattice. We have also performed calculations for particle 
numbers ranging from 30 to 80, for lattice sizes 8 3 and 10 3 , and various temperatures: 
in all cases, the results agree to within the aforementioned errors. 

According to the theory l53l l54l the asymptotic behavior in the limit of large 
momenta n(k) °< C(kf/k) n should at all temperatures be governed by the same 
exponent, namely n = 4. This is consistent with a value of the exponent n = 4.5(5) 



extracted from the MC data. Both the energy 9.1 and the entropy 9.2 exhibit a definite 
transition between low and a high temperature regimes separated by a characteristic 
temperature To: 

7b = 0.23(2)e f . (9.48) 
We shall discuss the relation between 7b, the superfluid critical temperature T c , and 



the pair breaking temperature T* in Sec. 9.2.8 First we focus on the low temperature 
limit. 

At T = 0, several interesting quantities describe the symmetric unitary system: one 
is the energy as expressed through the Bertsch parameter | = E$f I Efq; related is 
the somewhat fictitious energy of the interacting normal state §y — En /Efq; finally, 
there is the pairing gap A = tjef- The T = value of these quantities have been 
obtained to high precision by other groups using the variational fixed-node Monte 
Carlo techniques lfT8ll531l56ll57l . Unlike our approach, these T = techniques suffer 
from a sign problem that is overcome by using a fixed-node constraint: This formally 



provides only an upper bound on the energy. Our result £ = 0.37(5) (see Table 9. 1 
agrees with these variational bounds, B, = 0.44(1) HUE), % = 0.42(1) ll56l [57T . 
and with more recently quoted AFQMC results % = 0.40(1) |58l|59]. Although not as 
precise, our method is truly ab initio and hence provide a non-trivial validation of 
these variational results. 

The quantity for the normal state — though not precisely defined (since the 
normal state is not the ground state) — provides a useful description of the physics. 
For example, in the high temperature regime T > To, the energy is described well 
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by the energy of a free Fermi gas shifted down by 1 — <jj# (shown as a solid line 



in Fig. 9.1 1, where §y = t, + 8^ ~ 0.52 can be found by determining what shift is 
necessary to make the solid curve coincide with the high temperature data (where 
the gas is expected to become normal). 

Taking E, « 0.4 this gives the condensation energy 8^ « 0.12 which is roughly 
consistent with the estimate 



8E 5 ( A \ 2 



-0.15 (9.49) 



* fe F N 8\e F J 

based on the BCS expression for SE = %jpN (see [60]) and the QMC value of 
the pairing gap where A ~ 0.50£f ifTSl |57 | and confirmed by us in ISTIl (which 
turns out to be very close to the weak-coupling prediction of Gorkov and Melik- 
Barkhudarov ifTSl l62l ). Our estimate should also be compared with the results 
<^v ~ 0.54 of lfT8ll63ll and 4jv ~ 0-56 of [64 1 obtained by considering only normal 
state nodal constraints. Finally, a similar result §jy ~ 0.57(2) (see Eq. ( 9.95e] >) arises 



from fitting the SLDA density functional to be discussed in Sec. 9.3.2.1 

At low temperatures, T < 7b, temperature dependence of the energy can be ac- 
counted for by the elementary excitations present in the superfiuid phase: boson-like 
Bogoliubov- Anderson phonons and fermion-like gapped Bogoliubov quasiparticles. 
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Fig. 9.2 The entropy per particle with circles for 8 3 lattice, and with a dashed line the entropy of a 
free Fermi gas with a slight vertical offset. The statistical errors are the size of the symbol or smaller. 
From (45). 



Their contributions are given by 



£ph+qp(n = -E F N 




27r4 3 r 

— ^exp 



, (9.50) 
(9.51) 



The sum of the contributions from these excitations is plotted in Fig. 9.1 as a dashed 
line: Both of these contributions are comparable in magnitude over most of the 
temperature interval (To/2, 7b). Since the above expressions are only approximate 
for T <C T c , the agreement with our numerical results may be coincidental. 

At T > T c the system is expected to become normal. If 7b and T c are identified, 
then the fact that the specific heat is essentially that of a normal Fermi liquid Ep(T) 
above 7b is somewhat of a surprise: one would expect the presence of a large fraction 
of non-condensed but unbroken pairs. Indeed, the pair-breaking temperature has 
been estimated to be T* ~ 0.55£p, based on fluctuations around the mean-field, 
see H] 12 |65l [3] U |66]|. This implies that for T c < T < T* there should be a 
noticeable fraction of non-condensed pairs. In the next sections we will show that this 
is indeed the case and that above the superfluid critical temperature, the fermionic 
spectrum still contains a gap, giving rise to the so-called pseudogap phase. 

From the data for the energy E and chemical potential fi, one can compute the 
entropy S using the unitary relation PV = %E (true of a free gas as well) which holds, 
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where P is the pressure, V is the volume and E is the energy. It is straightforward to 
show that 

S _ E+PV-/XN _ %(x)-£(x) 
N ~ NT ~~ x ' 
where = h/Ef and x = T /£p determines the entropy per particle in terms of 



(9.52) 



quantities extracted from our simulation. As shown in Fig. 9.2 the entropy also 
departs from the free gas behavior below Tq. 

This data can be used to calibrate the temperature scale at unitarity |fl9l ISTl . 
Indeed, extending the suggestion of l67l . from a known temperature in the BCS 
limit, the corresponding S^T^cs) can be determined. Then, by adiabatically tuning 
the system to the unitary regime, one can uses S(T BC s) = S(7unitary) to determine T 
at unitarity. (In practice the experimental procedure goes in the opposite direction, 
namely measurements are performed at unitarity, and then the system is tuned to the 
deep BCS side, see ll68l .) 

On the other hand, knowledge of the chemical potential as a function of temper- 
ature allows for the construction of density profiles by using of the Local Density 
Approximation (LDA) (see the next section). In turn, this makes it possible to deter- 
mine S(E) for the system in a trap, fixing the temperature scale via dS/dE = l/T. 
Direct comparison with experiment shows remarkable agreement with our data (we 



discuss this later in Fig. 9.6 1 




Fig. 9.3 The critical temperature T c (squares either error bars) and the characteristic temperature 7o 
(circles with error bars) around the unitary point determined in QMC and using finite size analysis. 
On the far left BCS side of the critical point we show (solid green line) the expected BCS critical 
temperature, including the corrections due to induced interactions 11311621 . and on the far right side 
of the BEC side of the unitary point we show (solid green line) the expect critical temperature in the 
BEC limit. For more details see |45 1. 
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In the following we present a brief summary of our results near unitarity on both 



the BCS a < and BEC a > sides, (see Fig. 9.3 and Table |9.1| . The coupling strength 
was varied in the range —0.5 < l/kpa < 0.2 (where kp — (37t 2 n) '/ 3 ), limited on the 
negative (BCS) side by the finite volume V (which becomes comparable to the size of 
the Cooper pairs), and on the positive (BEC) side by the finite lattice spacing / (which 
becomes comparable to the size a — O(Z) of the localized dimers, manifesting as 
poor convergence of observables). 



1/kfU 


E(0)/E F 


To 




Eo/E F 


T c < n c /e F E c /E p 


-0.5 


0.60(4) 


0.14(1) 


0.685(5) 


0.77(2) 




-0.4 


0.59(4) 


0.15(1) 


0.65(1) 


0.75(1) 




-0.3 


0.55(4) 


0.165(10) 


0.615(10) 


0.735(10) 


0.105(10) 0.61(1) 0.64(2) 


-0.2 


0.51(4) 


0.19(1) 


0.565(10) 


0.725(10) 


0.125(10) 0.56(1) 0.61(2) 


-0.1 


0.42(4) 


0.21(2) 


0.51(1) 


0.71(2) 


0.135(10) 0.50(1) 0.54(2) 





0.37(5) 


0.23(2) 


0.42(2) 


0.68(5) 


0.15(1) 0.43(1) 0.45(1) 


0.1 


0.24(8) 


0.26(3) 


0.34(1) 


0.56(8) 


0.17(1) 0.35(1) 0.41(1) 


0.2 


0.06(8) 


0.26(3) 


0.22(1) 


0.39(8) 


0.19(1) 0.21(1) 0.25(1) 



Table 9.1 Results for the ground state energy, the characteristic temperature Tq, and the corre- 
sponding chemical potential and energy, from the caloric curves E(T) and the upper bounds on 
the critical temperature T c from finite size scaling and the corresponding chemical potentials and 
energies 1451 . 



9.2. 7 Response to External Probes and the Spectral Function 

In order to get an insight into basic degrees of freedom which contribute to the low 
energy excitations of the system one has to investigate the response of the system to 
various external probes. Here we will present the simplest possible probe: adding 
a particle to the system and calculating the probability amplitude of finding it in a 
given single particle state. This requires calculating the one-body finite temperature 
(Matsubara) Green's function (701 : 



Sf (p, x)= X - Tr{exp[-(/3 -t)(H- hN)] y f (l>) exp[-T(H - /itf) V(p)]}> (9.53) 

where j3 = l/T is the inverse temperature and T > 0. The trace is performed over the 
Fock space, and Z = Tr{exp[— j8 (H — /J.N)]}. The spectral weight function A(p, co) 
can be extracted from the finite temperature Green's function using the relation: 

^(P,t) =-~ T dcoA(p,co) i ex P(- CTT ) (9.54) 
v ' 2kJ-oo y * y l+exp(-coj3) 

By definition, A(p, co) fulfills the following constraints: 
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A(p,<o)>0, £^A(p,a>) = l. (9.55) 

Since our study focuses on the symmetric (unpolarized) system and the Hamilto- 
nian is symmetric under a -f-> b, £f (p, t) is block diagonal and the species index is 
suppressed in all formulae. The numerical evaluation of the one-body temperature 
propagator ( |9.53[ l is performed as described above, using a Trotter expansion of 
exp[— x(H — fiN)] followed by a Hubbard-Stratonovich transformation and Metropo- 
lis importance sampling. Details can be found in [61 ]. 



The numerical determination of A(p, co) by inverting ( 9.54 1 is an ill-posed problem 
that requires special methods. We have used two, based on completely different 
approaches. The first approach is the maximum entropy method 17T1 l72l l73l [741 . 
which is based on Bayes' theorem. Quantum Monte Carlo calculations provide us 
with a discrete set of values S?(p, T/), where i = 1,2,..., JV X = 50. We treat them as 
normally distributed random numbers around the true values Sf (p, T,). The Bayesian 
strategy consists in maximizing the posterior probability 

P(A\G) °c P(G\A)P(A) (9.56) 

of finding the right A(p, co) under the condition that £f (p, T/) are known. Here, 

P(G\A)~exp(-l-x 2 ) (9.57) 



2' 



is the likelihood function, where 



^ 9 

Z 2 = E [&(j>,Ti)-&(p,*i)] /o 2 . (9.58) 
The quantity ^(p, T,-) is determined by the spectral weight function in the discretized 



form of (9.54 1 at frequencies Cty. The prior probability P(A), describing our ignorance 
about the spectral weight function, is defined as P(A) <x exp(a5(.#)), where a > 
and S{^f) is the relative information entropy with respect to the assumed model 



k 



A(p,aJst)-^(03!fc)-A(p,fijfe)ln 



A(p,a>*A 

M (co k ) J 



(9.59) 



Hence the maximization of P(A\G) leads in practice to the minimization of the 
quantity l% 2 — aS{J() with respect to A lISTI . 



The second approach is based on the SVD of the integral kernel Jtf of (9.54 1, 
which can be rewritten in operator form as 

Sf(p,TO = (JfA)(p,Ti). (9.60) 

The operator possesses a singular subspace 

Jtifui = J^*\i = XiUj, (9.61) 
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where J€* denotes the adjoint of ', A, are the singular values, and m ; - and v, are 
right-singular functions and left-singular vectors respectively. The singular subspace 
forms a suitable basis for the expansion of the spectral weight function 1175117611771 . 
which we can then write as 

A(p,a)) = X>(pH(®), MP) = f (f(p)-Vi). (9-62) 

where (_ ■ _) is a scalar product and r is the rank of the operator Jfc$f* . Since (p, T,-) 
is affected by Monte Carlo errors a,-, the coefficients carry some uncertainty 
Abj. Each set of expansion coefficients £ (fc; — Abi,bj + Abj) reproduces Sf (p, T,-) 
within its error bars. We use this flexibility of choosing the expansion coefficients to 
produce a solution satisfying constraints (9.55 I |78 1. The relative advantages of each 



method will be discussed elsewhere 1791 . 

A sample of calculated spectral weight functions at unitarity are shown in Fig. 9.4 
In order to characterize the quasiparticle excitation spectrum we have associated with 
the maximum of A(p, a>) the quasiparticle energy E(p): 




Eip) = ± ^{i^ +U -^) +A > (9 - 63) 

where the effective mass m*, the effective potential U, and the "pairing" gap A depend 



on temperature, and jx is an input parameter. In Fig. 9.5 we compare the spectrum 
of elementary fermionic excitations evaluated in [57], with the one extracted by us 
from our lowest temperature spectral weight function. 



9.2.8 The Pairing Gap, Pseudogap, and Critical Temperature 

In order to find the critical temperature for the superfluid-normal transition one has 
to perform the finite size analysis discussed in the previous section. Following this 
procedure, our data for the condensate fraction of the unitary Fermi gas indicates that 
T c < 0.15(1 )ef, considerably lower than the characteristic temperature Tq =0.23(2) 
found by studying the behavior of the energy and the chemical potential (see Fig.|9.3| 



and Table 9.1 1. Even though this result for T c is close to estimates by other groups 
(see e.g. [80 8T1 [82ll ). it should be pointed out that the experimental data of [68 1 
shows a distinctive feature in the energy versus entropy curve at a temperature close 
to T (see (69)). 

It is notable that both methods (the maximum entropy method and the SVD 
method) admit a "gapped" spectral function above the critical temperature T c : a 
situation commonly called a pseudogap. It characterizes the range of temperatures 
where the system exists in an exotic state which is neither normal, nor superfluid, 
and defies a conventional BCS description. Therefore the onset of pairing and su- 
perfluidity can occur at different temperatures. On the other hand, the pseudogap is 
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1 2 3 0.05 0.1 0.15 0.2 0.25 



e(p)/e F T/e p 

Fig. 9.5 Quantities extracted from the spectral weight function A(p, a>) at T = 0.l£f at unitarity 
(from 1611 V Left: Quasiparticle energies E(p) (squares). The line corresponds to the fit to |9.63} . 
The circles are the results of Carlson and Reddy 1571 . (See also Fig. |9.10| where the same data is 
used to fit the SLDA density functional.) Right: The single-particle parameters. One should note that 
while the effective mans and the self-energy show a very weak temperature dependence across the 
phase transition, the pairing gap halves in value at T c and vanishes around Tq, 



easy to understand in the BEC limit where stable dimers exist well above the critical 
temperature. This gives rise to a pseudogap phase, where the system share a BCS-like 
dispersion and a partially gapped density of states, but does not exhibit superfluidity. 
Several groups have been advocating various aspects of pseudogap physics in the 
unitary Fermi gas for the past few years [4, 66 . 83 84 85, 86ll . 

There have been several experimental attempts to extract the pairing gap in 
ultra-cold dilute Fermi gases l87l l88l |89l and a theoretical explanation of these 
spectra was given in 1901I9TI . It was later shown in [92, 93, 94 95 | that these initial 
interpretations of the rf- spectra ignored the strong final state interaction effects. 
Recent experimental measurement of pair condensation in momentum space and 
a measurement of the single-particle spectral function using an analog to photo- 
emission spectroscopy, directly probed the pseudogap phase and revealed its existence 
for 1 /{kfo) ss 0.15 [96|. Although this lies on the BEC side, there are indications 
that the pseudogap persists well into the unitary regime Il97ll98l . 

Our calculations show that the spectral function reveals the presence of a gap in 
the spectrum up to about T* as 0.20£f (see Fig. 9.5 1, and a two peak structure around 
the Fermi level at temperatures above T c 0611 1991 . We note that T* is close to Tq 
(not surprising in hindsight), the temperature at which the caloric curve E(T) has a 
shoulder OH ED (called 7b in |05]). 
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9.2.9 Describing Trapped Systems with Quantum Monte Carlo 
Results 

The Monte Carlo calculations presented above assume that the system is uniform. 
In experiment, however, this condition is not fulfilled since atoms are trapped in an 
external potential which induces inhomogeneity of density distribution. Most of the 
atomic trapping potentials used in these experiments can be approximated rather 
well with harmonic potential wells. Such potentials can be shown to satisfy the virial 
theorem at unitarity, namely E(T,N) = 2N(U) — 3m<»?(z 2 ) |100|, and therefore 
simply measuring the spatial shape of the cloud allows for a unique determination 
of the unitary gas energy at any temperature. One of the main goals is therefore 
to provide a link between the results of experiment [68| and the available finite 
temperature QMC calculations. 

At unitarity (1/kpa = 0) the pressure of a homogeneous unitary gas is determined 
by a universal convex function hr(z): 

where T and jj, are the temperature and the chemical potential, respectively. 3^(T,n) 
is a convex function of its arguments (second law of thermodynamics) if and only if 
hj{z) is convex. One can show 1 1 1 1 that thermodynamic stability implies positivity 
hr{z) > 0, monotonicity h' T (z) > 0, and convexity hj(z) > 0. Remembering that 
the grand canonical potential is Q(V, T,fi) = —V3P(T,\l) one can show that the 
energy of the system reads: E = 33 g V /2, where V is the volume of the system. As 
it was mentioned before, this relation between energy and pressure is identical in 
form to the one corresponding to non-interacting particles. In the high-temperature 
limit ji — > — °o and £P(T,jj.) tends from above to the free Fermi gas pressure. In the 

low-temperature limit £?(T,}i) tends from above to &(0,E,Ef) — 4j5e 5 / 2 ^ /5. 

Standard manipulations show that all the thermodynamic potentials for the unitary 
Fermi gas can be expressed in terms of a single function of one variable, a property 
known as universality |[T9ll2Tll5Tll80l[8T1 . This property was incorporated in our 
interpolation. At high temperatures we notice that our results smoothly approach 
the corresponding free Fermi gas results with some offsets for the energy, chemical 
potential and entropy |[T9l |5D . 

At this point we assume that the Local Density Approximation (LDA) can be 
used to describe the properties of an atomic cloud in a trap. We will neglect the 
gradient corrections as one can show that for the mostly-harmonic traps used in 
typical experiments the role of the gradient corrections is relatively small [69], as the 
average interparticle distance, and thus the Fermi wave length, is much smaller than 
the harmonic oscillator length. 

In this approach, the grand canonical thermodynamic potential for a unitary Fermi 
gas confined by an external potential U (r) is a functional of the local density n(r) 
given by 



Tin 



25 



9.2 - The Quantum Monte Carlo Approach 



Q. = d'r 



where 



e F (r)p(*)n(r) +C7(r)n(r) - An(r) 



ft 2 



, ef (r) = -^[37r 2 «(r)] 2 / 3 , 
£/r (r) 2m 



(9.65) 



(9.66) 



and we have used the universal form for the free energy per particle F/N in the 
unitary regime: 



(9.67) 



where for a homogeneous system ^ (x) = 5E/3e F N, a(x) = S/N is the entropy per 
particle and x = T /e F . The overall chemical potential X and the temperature T are 
constant throughout the system. The density profile will depend on the shape of the 
trap as dictated by 8Q / 8n(r) = 0, which results in: 



8Q _8(F-XN) 



8n(r) 8n(r) 



= /i(x(r)) + t/(r)-A = 0. 



(9.68) 



At a given T and X, ( |9.66 1 and (9.68 1 completely determine the density profile n(r) 
(and consequently both E(T,N) and S(T,N)) in a given trap for a given total particle 
number. The only experimental input we have used is the particle number, the trapping 
potential and the scattering length at B — 1200 G, taken from |68]. The potential was 
assumed to be an 'isotropic' Gaussian, as suggested by the experimental group [68 1, 
although it is not entirely clear to us to what extent this is accurate, especially in the 
axial direction. We have approximated the properties of the atomic cloud at B — 840 G 
with those at unitarity (B = 834 G), where we have MC data. For B = 840 G and for 
the parameters of the Duke experiment [68 1 one obtains 1 jk F a = —0.06, using data 
of H102L if the Fermi momentum corresponds to the central density of the cloud at 
7 = 0. 

Our results for the entropy of the cloud and the density profiles for several 



temperatures, are shown in Figs. |9.6| and 9.7 In all the figures the temperature is 
measured in natural units of £ F (0), corresponding to the actual central density of 
the cloud at that specific temperature. In J68l 11031 the temperature is expressed 
in units of the Fermi energy at T — in a harmonic trap: e F " = hQ(3N) 1 ' 3 . It is 



clear from Fig. 9.7 that the central density decreases with T and that the superfluid 



core disappears at T c = 0.23(2)£f(0), which translates into T c = 0.27 (3)£p to be 
compared to T c = 0.29(2)4° of EQ. There is a noticeable systematic difference 
between theory and experiment at high energies, see Fig. |9.6| This discrepancy can 
be attributed to the fact that the experiment was performed slightly off resonance, on 
the BCS side, where l/k F a = -0.06. 

Recently a couple of new experiments have been published, one by the Paris 
group [104] and another by the Tokyo group [105|. Using new techniques these 
groups were able to extract directly from cloud images the pressure as a function of 
the fugacity. While the Paris group has observed a very good agreement with our 
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(E(T)-E(0))/E 

Fig. 9.6 Entropy as a function of energy for the unitary Fermi gas in the Duke trap 1 68 1 : experiment 
(points with error bars) and present work (solid curve), where E$ = N£p°. Inset: log-log plot of 
E(T) as results from our calculations and as derived from experimental data |68|. The temperature 
is units of the corresponding Fermi energy at the center of the trap: £f(0). From |69|. 




Fig. 9.7 The radial (along shortest axis) density profiles of the Duke cloud at various temperatures, 
as determined theoretically using the QMC results 1 19 ST). The dotted blue line shows the superfluid 
part of the cloud, for which x(r) = T/ef(i) < 0.23. The solid red line shows the part of the system 
that is locally normal. Here aj m = h/mCO max . From l? 69l . 
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Fig. 9.8 The comparison between the ration of the pressure versus fugacity of a unitary Fermi gas 
and the pressure of a free Fermi gas: measured in 1 106 1 (red filled circles) and calculated in QMC 
in 1 69 1 (blue filled diamonds). 



QMC results, see Fig. 9.8 they have also noticed that the results of the Tokyo group 
show systematic differences [106|. 

One can summarize that so far the bulk of the theoretical predictions obtained 
in ab initio QMC have been confirmed experimentally with impressive accuracy in 
most cases, often at a level of a few percent, which is the accuracy of both theoretical 
calculations and of many experimental results as well. The emergence of a pseudogap 
in the unitary gas is a fascinating new feature, but still in its infancy both theoretically 
and experimentally. 



9.3 Density Functional Theory for the Unitary Fermi Gas 

The idea of Density Functional Theory (DFT) originated with Hohenberg and 
Kohn 1110711 and Kohn and Sham lfT08l (see the monographs lflQ9l [TTOl for an 
overview) where they proved that the ground state energy and the density of a system 
of interacting fermions in an arbitrary external potential V e xt( r ) ma Y be found by 
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minimizing a functional 

£[n(r)]+|d 3 rV ext (r)«(r). (9.69) 

The utility of this approach is that the functional E[n(r)] depends only on the inter- 
actions of the system and is independent of the external potential. Thus, if we were 
able to deduce E[n(r)} for the unitary Fermi gas, then by simply minimizing a single 
functional, we could determine the ground state in any external potential, including 
arbitrary trapping geometries and optical lattices. 

The challenge is that the Hohenberg-Kohn theorem is an existence theorem. The 
exact form of the functional E[n(r)] is unknown, and in general it may be extremely 
complicated and highly non-local. In problems that are under perturb ative control, the 
functional can be formally derived (see [ 1 1 1 J), but in highly non-perturbative prob- 
lems such as the unitary gas, one must choose a physically motivated approximate 
functional and check its accuracy. 

Our strategy is thus: 

1. Postulate simple functional forms capturing the relevant physics with a small 
number of parameters. 

2. Use ab initio results to fix these parameters. 

3. Validate the functional with different ab initio and experimental results. 

4. Make interesting and verifiable physical predictions. 

The computational cost of minimizing the density functional is much less than 
solving for many-body wavefunctions, and one may consider substantially larger 
systems, untenable with ab initio methods. This allows one to make direct contact 
with typical mesoscopic experiments for example. In this way, one may view the 
density functional as a bridge between microscopic and mesoscopic physics. 

As we have noted, although DFT is exact in principle, for non-trivial systems 
we must postulate a form for the functional. Nevertheless, it provides a substantial 
improvement to the ad hoc mean-field methods typically employed to study the 
properties of large non-perturbative many-body systems. Without a program for 
systematically correcting the functional, the DFT approach will not be the final word. 
However, judging from the success of the approach in quantum chemistry, and from 
the results presented here, we expect that without too much effort one should be 
able to obtain percent level accuracy for a wide range of systems, which should be 
sufficient for quite some time. 

The qualitative success of the Eagles-Leggett [ 1 , 2 1 mean-field model describing 
the BCS-BEC crossover suggests a functional description of the unitary Fermi gas 
in terms of quasi-particle fermionic states (see ( |9.76[ >). As discussed in Sec. |9.1| 
although the BdG approximation is quite successful, it is quantitatively inaccurate as 
it describes all interaction effects through the condensation energy (pairing) alone, 
completely omitting the "Hartree-Fock" contribution which dominates the energetics. 
To see this, consider the typical local interaction gd'b^ba between species a (spin 
up) and species b (spin down). The mean-field approximation retains the pairing 
term g (a'b^) (ba) = gv^v and the Hartree term g (a* a) (b^b) = gn a rib. (The other 
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quadratic Fock term (a T b) (b^a) has zero expectation.) The problem arises upon 
renormalization: As discussed below, the anomalous density v is formally divergent, 
and regularization requires taking the coupling g —> to keep the gap parameter 
A = —gv finite. Since the densities remain finite, the Hartree contribution gn a rib — > 
vanishes. 

In weak coupling, one can carefully take the zero-range limit while summing 
ladders lfl2l ITOll . obtaining the well known form an a nb of the Hartree interaction, 
which is clearly invalid in the unitary limit \a\ — > °°. In particular, for the symmetric 
phase n a = rij, = n, there is no additional length scale, and so we must have a 
dependence ~ n 5 / 3 as dictated by dimensional analysis. This physics — the dominant 
contribution to the energetics (see the discussion below ( |9.49[ l) — is completely 
missing from the BdG (mean-field) approach and is one of the main deficiencies we 
hope to overcome within an improved DFT description. 

We shall first discuss an improved local DFT for symmetric systems n a = tif, = n '- 
the Superfluid Local Density Approximation (SLDA). This is a generalization of 
the Kohn-Sham Local Density Approximation (LDA) to includes pairing effects and 
subsumes the BdG form, adding an n 5 / 3 Hartree interaction term. 

We subsequently extent the SLDA to study asymmetric systems n a ^ rib through 
the use of the Asymmetric SLDA (ASLDA) functional that subsumes the SLDA. 
The approach of both these approximations is to introduce as few parameters as 
possible that are consistent with the scaling and symmetries of the problem, then 
to determine the coefficients of these terms by matching to ab initio properties in 
the thermodynamic limit. The form of the functionals is described in Sec. 9.3. 1 the 
fitting of the parameters is discussed in Sec. 9.3.2 and some physical applications 



are presented in Sec. 9.3.3 



9.3.1 The Energy Density Functional 

We start with the most restrictive conditions of a cold (T = 0) symmetric (n a = rij,, 



m a = nib) unitary (|a| = °°) Fermi gas. As discussed in Sec. 9.1 the only dimensionful 
scale in the problem is the density n, so dimensional analysis provides significant 
constraints on the form of the functional and thermodynamic functions, allowing 
us to postulate a simple functional form characterized by only three dimensionless 
parameters. Relaxing any of these conditions will introduce additional dimensionless 
parameters. In particular, we consider the dimensionless polarization p — (n a — 
rib)/(n a +rib) to formulate ASLDA II101II1121I . The generalized ASLDA functional 
promotes the dimensionless parameters to dimensionless functions of this asymmetry 
parameter p. 
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9.3.1.1 Local Density Approximation (LDA) 

In general, the energy functional might be a highly non-local and extremely com- 
plicated object. One major simplification is to assume that the functional is local 
and can be represented by a function of various types of densities. This amounts to 
introducing the energy density S which is a function (as opposed to a functional) of 
the local densities and their derivatives (referred to as gradient corrections): 

E KS = |d 3 r^ s [«(r),T(r),V«(r),...]+f/(r)«(r) + ---, (9.70) 

where U (r) represents an external (trapping) potential. This local density approxi- 
mation (LDA) has met with remarkable success in quantum chemistry applications 

H113lll09lfTT0l . 

The simplest function contains a single term E °< n 5 / 3 . This — along with gradient 
corrections — has been explored in 11 1 141 II 151 . and, while it can model the energetics 
of the symmetric gas, it does not include information about pairing correlations. The 
extensions we describe here include both kinetic terms and an anomalous pairing 
density. 



9.3.1.2 Densities and Currents 



The first task is to construct the densities and currents. In the SLDA, we con- 
sider five types of densities: the standard particle densities n a {r) ^ (a^(r)a(r)) 
and «i(r) °< (b^(r)b(r)), the kinetic densities T fl (r) °< (a^ (r)Aa(r)) and T^(r) °c 
(b' (r)Ab(r)), an d an anomalous density v(r) °< (a(r)Z?(r)). When considering time 
dependence (Sec. |9.4|l, we must also include the currents j„ (r) °< (a T (r)V a(r)) a nd 



9.4.2 



In 



jf>( r ) 06 (^(r)Vb(r)) to restore Galilean invariance as discussed in Sec 
principle, these densities may be non-local, but to simplify the functional we wish to 
consider only local quantities. The local form of the anomalous density v leads to 



UV divergences that we must regularize as we discuss in Sec. 9.3.1.4 



The formal analysis proceeds with a four-component formalism discussed in 



Sec. 9.6.2 but the symmetries of the cold atom systems allow everything to be 
expressed in terms of two-component wavefunctions (see Appendix|9.6|l 



Wn{r) 



u„(r) 
v„(r) 



(9.71) 



with energy E n . The densities and currents are constructed from these as 
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» a (r) =£K(r)| 2 //3(£«), n b (r)=^\v n (r)\ 2 M-E n ), 

" " (9.72a) 

T„(r) =£|V«„(r)| 2 //3(£„), T,(r) =^|Vv„(r)| 2 //3(-£ n ), 

V(r) = ^Z U n( r K( r )(fp(- E n)-fp(En)), 

jflW = ^£KMVw n (r)-%(r)VM*(r)]/0(£„)^ (9.72b) 
J»W = ^I[v:(r)Vv„(r)-v n (r)Vv:(r)]/^(-£„), 

Z n 

where fp(E n ) = l/(exp(j3£„) + 1) is the Fermi distribution and j3 = l/T is the 
inverse temperature. Even though we shall only discuss the zero temperature limit of 
SLDA it is convenient for numerical purposes to introduce a very small temperature 
(much smaller than any other energy scale in the system) so that S'(li) is a smooth 
function. 




9.3.1.3 Functional Form 

Our functionals generically include a kinetic term and a pairing term of the form 

-^+gv f v + ---, (9.73) 

along with additional density dependent terms, where all of the densities and cur- 
rents n(r) etc. are functions of position but have no non-local structure. (Note that 
here and in many of the following formulae we suppress the explicit dependence 
on position r.) In the superfluid, this local approximation has formal difficulties 
since the anomalous density v(r,r') ~ £M„(r)v*(r') ~ |r — r'|~' diverges for small 
|r — r'| if the pairing field is taken to be a multiplicative operator A (r). The kinetic 
energy densities T a .b( r ) diverge as well. A proper local formulation thus requires 
regularization 11 161 as discussed in Sec. 9.3.1.4 We introduce an energy cutoff 



E c — V c (r) ~ L|£|<£ c M n( r ) v n( r ) — an d a cutoff dependent effective interaction g e g 
such that 

A = -gv = -geffVc (9.74) 

is finite and independent of the cutoff as E c —> °°. Once this is done, we can write the 
functional as 

T + lb 




-A'v + ---. (9.75) 

Note that v is still formally divergent, but will cancel with a similar divergence in 
the kinetic piece such that the energy density is finite. The full forms of the local 
functionals considered here are thus: 
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Bogoliubov de-Gennes (BdG) lfIT7l : 

^G^P^ + P^+gv'v. (9.76) 
2m a 2m b 

For homogeneous systems, this is equivalent to the Eagles-Leggett mean-field 
theory where the parameters here represent the bare parameters (elsewhere we 
shall only consider m a = m b = m) and the coupling constant is tuned to reproduce 
the vacuum two-body scattering length a. Note the absence of a self-energy: 
all of the interaction effects are modelled through the pairing interaction. One 
unphysical consequence is that the normal state is described as completely non- 
interacting in this model. While this may capture some qualitative features of 
the theory, and provides a rigorous variational bound on the energy, it cannot be 
trusted for quantitative results beyond the rather poor variational upper bound. 
SLDA: 



h 2 fa, N „ 3 



'SLDA 



(T„ + T i )+j3^-(3^ 2 ) 2 / 3 K+»,) 5 / 3 +gyV (9.77) 



w \ 2 r 10 



This may be thought of as the unitary generalization of the symmetric BdG 
functional for symmetric matter n a = n b = n + /2 to include a self-energy term 
n 5 / 3 (whose form is fixed by simple dimensional analysis) and an effective mass 
'«eff — m/a. The three parameters here a, j3, and the pairing interaction g must 
be fixed by matching to experiments or ab initio calculations as discussed in 



Sec. 9.3.2.1 Since g is formally zero in the large coupling limit, we characterize 
it with a dimensionless constant y such that g^ — (n a + n b ) 1 / 3 /j— A where A is 
the cutoff discussed in Sec. 19.3. l~4l 
ASLDA: 

^aslda = — (a a (n a ,n b )^- + ab(n a ,n b )—+D(n a ,n b )) +gV f V. (9.78) 

Here we allow for polarization n a ^ n b and so we must generalize the param- 
eters such as the effective masses and self-interaction to be functions of the 
local polarization p = (n a — n b )/(n a +ni,). Dimensional analysis restricts these 

a a ,b{kn a , Xni,) = a(n a ,n b ) and D{Xn a ,Xn b ) = X 5 ^D(n a ,n b ) so that we need 



only to parametrize functions of the single variable p as discussed in Sec. 9.3.2.1 



To fully define these functionals, we must now regularize the pairing interaction g 



(Sec. 9.3.1.4 1 and then specify the values and functional forms of the parameters and 



parametric functions (Sec. |9.3.2| . 



9.3.1.4 Regularization 

As formulated, the local theory is ultraviolet divergent due to the well known be- 
haviour of the anomalous density: 
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v(r,rV£ Mn (r)v*(r')°c^— r (9.79) 

There are many ways of dealing with this. For example, physical potentials are always 
non-local, and the non-locality naturally regulates the theory. However, in the unitary 
gas, the non-local (range of the interaction) is much smaller than any other length 



scale in the system and the stability of the system (see Sec. 9.1 1 indicates that the 
low-energy large-distance physics should be independent of the short-range details. 

As a result, one can choose any sort of regularization scheme that is convenient 
and obtain the universal results with an appropriate limiting procedure. In the homo- 
geneous case, one can use a variety of techniques: some interesting choices include 
dimensional regularization 111 181 and selective distribution functions [ 1 19 1. The most 
straightforward is to use a momentum cutoff, but for inhomogeneous systems, mo- 
mentum is not a good quantum number. Instead, an energy cutoff E c suffices. All 
quantities — especially the divergent anomalous density — can be computed from 
states with energies below this cutoff: 

v- *fp( E n)-fp(-En) 

v c = L u " v n if- ■ ( 9 -80) 

(To improve the behaviour, we actually use a smooth cutoff so that discontinuities 
are not introduced when levels cross in and out of the sum during the self-consistent 
iterations.) 

To better understand the nature of these divergences, consider the ultraviolet limit 
where the length scale is much smaller than any other scale in the system. In this 
limit, the semi-classical Thomas-Fermi approximation may be applied locally. The 
linear divergences in both the symmetric combination of the kinetic energy and in 
the anomalous density have the form 

2(m*) 2 A f A m*A 

T + {k) = T a {k) + T b {k)^ Z{m ^ 2 V ^^^ (9 ' 81) 

where the average effective mass m* =m/a+ = 2m /{a a + 0%) enters explicitly 
through the equations of motion. From this it is clear that the combination 

h 2 x+ t h 2 (a a x a a h T b \ t 
2^- A V =m(— + — ) +8V V 

remains finite if we regularize the theory such that the gap parameter remains finite 
for all values of the cutoff 

A = - 8eB V c . (9.82) 



When regularizing the BdG equations (9.76 1, we hold fixed the vacuum two-body 
scattering length, 

m _ 1 1 f d 3 k 1 

47z¥a ~g + 2j j2nT¥k 2 7 (9 ' 83) 
— — +i0+ 
Zm 
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where j is the principal value integral. This may be easily derived from the pseudo- 
potential approach (see for example ||120||12ll or for higher partial waves B122I 1 123D , 
In the other DFTs \9.11\ and ( |9.78| l, g does not represent the physical interaction, 
but is simply another parameter of the theory. Thus, we define a similar regulariza- 
tion scheme by introducing a finite function C(n a ,np) that must be fit in order to 
characterize the pairing interaction and correlations^ 

a+v 1 f d 3 k 1 a+ 

C(n a ,n b ) = -^ + -j— 3-2-2 = ^+A. (9.84) 

A 2 J {2%) ^l_M+ +i0+ £eff 

2m a+ 



This differs from (9.83 i in two ways: 1) we have included a factor of the effective 



mass parameter a+ to ensure that the divergences (9.81 1 cancel and, 2) we have 
shifted the pole of the integral by the average local chemical potential fi + = (jj, a — 
V fl + Hb+Vb)/2 to improve convergence. As pointed out in 11 161 . the shift does not 
change the integral in the limit of infinite cutoff, but greatly improves the convergence 
if a cutoff is used. Given a fixed momentum cutoff k <k c , the integral A in the second 
term can be performed exactly 

m k c f k k c + k \ 
A = ¥^{ 1 -2k- ln k-^- \ (9 ' 85) 

where h 2 ^/ (2m) = jx + /a+ defines the location of the pole. In general, translational 
invariance is not preserved, and so we must use the fixed energy cutoff \E(k) \ < E c 
that enters ( |9.80| l rather than a momentum cutoff as the latter is not a good quantum 
number. To relate the two we used the local quasiparticle dispersion relationship: 

h 2 

— a+(r)^( r )- M+ (r)=0, (9.86a) 

^ n a + (r)k 2 (r)-fi + (r)=E c . (9.86b) 

This defines a position-dependent momentum cutoff k c (r) and effective coupling 
constant g(r) that can be used to regulate the anomalous density at any point in space: 



v ; h 2 2k 2 \ 2k c (r) k c (r)-ko(r) J 
° t+(rl -C(n a (r)Mr)) -A(r), (9.87b) 



5eff(r) 

A(r) = -^ eff (r)v c (r). (9.87c) 



1 We have changed notations slightly from 11241 using C(n a ,n^) = a+C(n a ,ni,) which simplifies 
the equations because, in the limit of infinite cutoff, A is independent of any densities and functional 
parameters. 
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Varying the functional with respect to the occupation numbers (see Appendix |9. 6. 1| 
for a formal description) allows us to derive the self-consistency conditions. Recall 
that the functional has the form 

tip" T ( T 2 ff \ /z^ 

a-(n a ,n h )—— + a+(n a ,n b ) — h— VjV c +— D(n a ,n b ), (9.88) 

2m \ 2m a+ J m 

and that, in the limit of infinite cutoff, A has no dependence on the functional 
parameters so thaj^] 

«+\ . / geff\ = _f 8eK\ 2 



dC = d — d — = - — dC. (9.89) 

VW \cc+J \a+J 

Thus, we have the following equations: 

Ka - Ha+Va A 1 " \ fu„\ = (u, 

A -K b + n b -V b \y n "U 



where 



h 2 

K a u = -—Vi(a a (n a ,n b )Viu) 
h 2 

K b v= - — V,-(a i (n a ,n t )V,-v) 



da,-(na,n b )h 2 %- da + {n a ,n b ) (h 2 X + 4 t V 



dC(n a ,n b ) A^A h 2 dD(n a ,n b 



OC±(n a ,n b ) = \[aa(n a ,n b )±a b (n a ,n b )], 

T± = T ±Tj,. 

and similarly with a-^rb. 



-Ua{r), 



9.3.2 Determining the 5LDA and ASLDA Energy Density 
Functionals 

Unless one has perturbative control over the theory, one cannot in general determine 
the correct functional from first principles. Instead, the functional must be treated 
as a model incorporating the most relevant physics for the application at hand. As 
such, one must determine some parameters in order to make predictions about other 



2 There is a small correction due to the residual density dependence of A at finite cutoff but in 
practice this is insignificant. 
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properties of the system. Here we use properties of homogeneous matter in the 
thermodynamic limit to determine the parameters of our functional, and then use 
the functional to compute the properties of non-uniform systems such as trapped 
gases. Our hope is that the single particle states in the self-consistent approach will 
provide a good description of the finite size (shell) effects missing in the Thomas 
Fermi approximation. 

Fortunately, the thermodynamic functions describing the unitary Fermi gas are 
tightly constrained 1 101 1, and have both calculational and experimental verification. 
We shall now describe how to use these constraints to determine the form of the 
dimensionless parameters describing the functional. 



9.3.2.1 Homogeneous Matter 

A simple Thomas-Fermi calculation can be employed to describe states of homoge- 
neous matter by exploiting the translational invariance of the system. This allows 
us to fix all non-gradient terms in the functional. The only remaining term — the 
effective mass — must be fixed by other means and we use the quasiparticle properties 
to determine this coefficient. 



Normal Phase 

The energy-density for the normal phase of homogeneous matter has the form 

5/3 



h 2 (6n 2 (n a +n b )) 



£\n a M = -- ' 2 ~ G(p), p = ^-^e[-l,l]. (9.90a) 



where 



G(p) = a( P )[ V /3 + a (-p) ( l —l) 5/3 +2 -y^ (p) (9.90b) 



2 

and P{p) is defined through 

5/3 



{bn 2 (n a + n b ) 



D(n a ,n b ) = ± 2Qn2 1 2- 2 '^{p). (9.90c) 

The function G(p) will be the main function that enters our numerical formulae]^] 

We shall define the dimensionless function G(p) by fitting a simple even polyno- 
mial to the Monte-Carlo data tabulated for /[/7(x)]nFrom this, the function D(n a ,n&) 



3 In our previous calculations 1 125 . 124], we used a more complicated parametrization: the present 
form G(p) is just as good and much simpler and we advocate its use instead. 

4 G(p) is related to the other dimensionless functions f(x) and g(x) discussed in the literature as: 
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may be directly expressed in terms of the inverse effective mass a{p), which may be 
independently parametrized: 



D(n a ,n b ) - 



(6K 2 {n a 



5/3 



207T 2 



G( P )-a(p) 



5/3 



-a(-p) 



5/3' 



The function G(p) describing the normal state has been well-constrained by 
Monte-Carlo data |64j (see Fig. |9.9| >. As shown in Fig. 9.9 the function G(p) is very 
well parametrized by a simple quadratic polynomial: 



G(p) =0.357 + 0.642//. 



(9.91) 



Symmetric Superfluid State n a = n b 

As suggested in 11271 . by considering the calculated properties of the fully paired 
symmetric superfluid, one may determine the values of the functions a(p), C(n a ,ni,), 
and D(n a ,n b ) at the point p = where the energy density functional depends only 
on the symmetric combination of parameters n a =n b and z a = z b . For any value of 
the inverse effective mass a = a„ = o, one can uniquely determine the self-energy 
/3 = j3 p= o and pairing interaction y by requiring that the energy and gap satisfy 

S SF = S{n,n) = 4<f f G = £ - ^y'J 3 , (9.92a) 

m 107T Z 

A = r,e F = ri- K ' . (9.92b) 
m 2 

The parameters 4 and tj = A /ep have been calculated using several Monte-Carlo 
techniques (l8l|57]|59]|6TJ|45|. We take the following values in our estimates Il57ll59l : 

- " { "'" ] -0.40(1), 77 = — =0.504(24). (9.93) 



S FG {n,n) e F 

In order to determine the effective mass, we consider the quasiparticle dispersion 
relationship l57l . Within our density functional, this has the form 



The function g(x) = g[p(x)\ introduced in 1 1 Oil has the necessary and sufficient requirement of 
convexity to satisfy the second law of thermodynamics; and the function f(x) = f[p(x)] was 
introduced in 11261 and has been tabulated using Monte-Carlo methods 1641 . 
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Fig. 9.9 Monte Carlo data used to fit the function G(p) (top) and in its raw form f(x) = g 5//2 (x) 
(bottom) representing the energy of the normal state with respect to the energy of the free system. 
We excluded the red points from our fit because we suspect that they slightly contaminated by the 
superfluid state (and hence have a lower energy). Fitting these close to the superfluid state would 
require a double hump structure in G(p) for which we do not yet see any physical motivation. To 
anchor the solution in the superfluid phase, we include a datum j8 p= o extracted from the symmetric 
state j9.95b| . The value here depends slightly on whether or not we also extract an effective mass, or 
hold a = 1 constant. Both fits are shown (but lie on top of each other). The present fit is the simple 
two-parameter quadratic given in {9.91} . At the lower-right of the lower plot we have shown the 
values of f x= \ for the superfluid state (black point). Finally, for comparison, we have included the 
function f(x) obtained using the standard mean-field (Eagles-Leggett) approximation as a dotted 
yellow line to show that it bears little resemblance to the physical curves. 
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£ - (&)= v(£" Meff ) +A1 (9M) 

where n a + rib = k F /3n 2 is the Fermi wave-vector and /z e ff is the effective chemical 
potential. It turns out that /i e ff also depends on A /Bp, so the quasiparticle dispersion 
relation is really a function of only two parameters: the effective mass and A/ef- 



The fit to the Carlson-Reddy data 012811 is shown in Fig. 9.10 and gives the 
following parameter values]^] 

a p=Q = m; f l/m- 1 = 1.094(17), (9.95a) 

j3 p=0 = -0.526(18), (9.95b) 

y" 1 = -0.0907(77), (9.95c) 

ri=A/e F = 0.493(12), (9.95d) 

E, N = a + J3 = 0.567(24). (9.95e) 

where ^n'S'fg lS the energy of the interacting normal state predicted by the functional. 



Note that this agrees very well with the value given by G(p) in Fig. 9.9 (we have 
used this parameter as an additional point in the fitting of G(p)). 

In principle, one should use some form of ab initio calculation or experimental 
measurement for polarized systems to determine the dependence of the parameters 
a, j3, and yon the polarization p = n\,jn u . Unfortunately, the fermion sign problem 
has made this difficult and there is presently insufficient quality data to perform such 
a fit. Instead, we simply fix 

Y(p) = 7p=o = -11.11(94). (9.96) 

If high quality data about polarized superfluid states become available, one might 
consider promoting this parameter to a polarization dependent function similarly to 
a(p) and G(p). This fixes the pairing interaction: 



m a + (p)(n a + n b y^ 

ti 2 rip) 



C(« fl ,n/,) = ^2 • ( 9 - 97 ) 



Effective Mass Parametrization: a(n a ,nb) 

As discussed above, the effective mass cannot be determined solely from the 
properties of homogeneous matter. It is also clear in DFT's developed perturba- 
tively 1 1291 1 130 1 that the effective mass is arbitrary. In the ASLDA, however, the only 
gradient terms that enter the functional are the kinetic terms T whose coefficients are 
the effective masses. To allay the need for additional gradient corrections, one must 



We have performed a simple two-parameter non-linear least-squares fit which has a reduced 
Xred = 1 ■ 1 ' indicating a very good fit. 
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provide a parametrization of the effective mass. Fortunately, three values are well 
determined: In a fully polarized system, the effective mass of the majority species re- 
mains unchanged, m p= \ — 1 .Qm, while in the minority species, the effective "polaron" 
mass m p= _ i = 1 .20m 11311 . We use this value, but note that there are other estimates: 
Monte Carlo calculations give m_i = 1.04(3) 164) and m_i = 1.09(2) 11321 , and 
experiments measure m_i = 1.06 (no error given) 11 1 3 3 H and m_i = 1.17(10) II 1 341 . 



The third value for symmetric matter mo = m/a p= o is determined in ( 9.95a I. 

We now have three data-points constraining the effective mass parametrization 
of a{p). For numerical reasons, in order to ensure that the effective potentials V a ^ 
approach zero as the density falls to zero, we impose the additional constraint that the 
first and second derivatives of a(p) vanish at the end-points p = ±1. Taken together, 
this fixes a sixth order, two parameter polynomial approximation for a{p): 

a{p) = 1.094 + 0.156p(l-2p 2 /3+p 4 /5)-0.532p 2 (l-/5 2 +p 4 /3). (9.98) 



9.3.2.2 Summary 

Here we summarize the complete definition of the ASLDA functional. The SLDA 
functional follows by setting the local polarization 

, . n„(r) — riiJr) 
P(r) = (- (9.99) 
n a {r)+n h (r) 
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l.l 



c 1.0 



II 

0.9 



0.8 




p = (n a -n b )/(n a +n b ) 



Fig. 9.11 Inverse effective mass a(p) = m/m e ff(p) as a function of the polarization p = (n a 
n b) I \ n a + i&) ■ The functional fit is the polynomial l|9.98|. 



to zero. First, fitting the quasiparticle dispersion relationships, gap and energy for the 
superfluid state gives the SLDA parameters at p = 0: 

a p=0 = 1.094(17), Pp=a = -0.526(18), y ; ;j = -0.0907(77). (from (gg5)) 

Using this derived effective mass, and the energy data for the normal state from 
Monte Carlo data we obtain the following polynomial fits defining the polarization 
dependence of the effective mass and self-energy: 



a(p) = 1.094 +0.156/7 ^1 - 



2p 2 - 4 



P_ 
5 



532/( 1-/ + ^ 



G(p) =0.357 + 0.642/? 2 , 

y(p) = yp=o = -ii.ii(94). 

These fix the specification of the functional parameters 



(from d9T98l ) 
(from fl£9l) ) 
(from ( [9^96) ) 



/ \ / \ , \ , . ~, . m a + (p)(n a + n b ) 1 / 3 

oc a {n a ,n b ) = a(p), a b {n a ,n b ) = a(-p), C{n a ,n b ) = — . 

fr 7{p) 



D(n a ,n b ) 



[6n 2 (n a 



5/3 



207T 2 

in terms of the densities 



G(p)-a(p) 



1+P 



5/3 



-a(-p) 



5/3' 
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n a (r)= £ \u n (r)\ 2 f p (E n ), n b (r) = £ \v„(r)\ 2 f p (~E n ), 

\E n \<E c \E n \<E C 

T fl (r) = £ \Vu n (r)\ 2 f p (E n ), r b (r) = £ \Vv n (r)\ 2 f p (-£„), 

|£„|<£ c |£„|<£ c 

Z |£„|<£ c 

j«W=' £ [u*(r)Vu n (r)-u n (r)Vu:(r)]fp(E n ), 

\E n \<E c 

m = \ £ K(r)Vv„(r)-v„(r)Vv:(r)]/ p (-£ n ), 

|£„|<£ c 

in the form 

^ a a (n a , n b ) y + a A (n a , n& ) y + D (n a , n b ) ) + g e ff V f V 



'ASLDA — 

m 



together with the renormalization conditions 



A(r) = -s eff (r)v c (r), = C(r) -A(r) 

m fe c (r) 

A(r) = 



/i 2 27T 2 



f k Q (r) k c (r)+k Q (r) \ 
\ 2k c (r) k c (r)-k Q (r)j- 



^ «+ (r)*g (r) - m+ (r) = 0, |- a+ (r)fc 2 (r) - M+ (r) = E c 

where a+ = (a a + a b )/2 and jti + = (jU a — V fl + ju^, + V b )/2 is the average chemical 
potential defined through the equations: 



where 



A -K b + fi b -V b I \v n \v, 



h 2 

K a u = - — Vi(a a (n a7 n b )Vju) 

h 2 

K b v = - — V i(a b (n ai n b )V iv) 

y _ da_{n ai n b )h 2 %_ | da + (n a ,n b ) f h 2 X + A^V 

dn a 2m dn a \ 2m a + (n a ,n b ) 

dC(n ai n b ) A^A h 2 dD(n ai n b ) 

■ + — 51 +Ua(T), 



dn„ a i m dn, 



a 



CC±(n a ,n b ) = 1 [a a (n a ,n b )±a b (n a ,n b )], 
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and similarly with a <-> b. 



9.3.3 Using the 5lda and Aslda 



Once the form of the DFT and its parameters have been fixed, the function needs to be 
tested and applied. Since we fit the parameters using QMC results for homogeneous 
matter, a non-trivial test is to compare it with ab initio results in inhomogeneous 
situations. This will asses the accuracy of the approximation we have made in 
neglecting gradient corrections beyond the kinetic term. In Sec. |9.3.3~T] we compare 
the predictions of the DFTs with QMC calculations of trapped systems. Next we show 
how the functionals can be used to explore mesoscopic physics inaccessible to QMC 
analysis techniques: we consider the structure of superfiuid vortices in Sec. |9. 3. 3. 2 



and the prediction of a supersolid phase in the asymmetric case in Sec. 9.3.3.3 




Fig. 9.12 The comparison between the GFMC 11351 , FN-DMC 11361 and SLDA total energies 
E(N). The clear odd-even staggering of the energies is due to the onset of the pairing correlations. 
The inset shows the discrepancy between the GFMC and FN-DMC and SLDA energies, SE(N) = 
Emc( n ) / e slda(N) - 1, where E MC {N) stands for the energies obtained in GFMC or FN-DMC 
respectively. 



44 



9.3.3 - Using the Slda and Aslda 



Normal State 



(N a ,Nb) EfNDMC 


E ASLDA 


(error) 


(3,1) 


6.6±0.01 


6.687 


1.3% 


(4,1) 


8.93±0.01 


8.962 


0.36% 


(5,1) 


12.1 ±0.1 


12.22 


0.97% 


(5,2) 


13.3±0.1 


13.54 


1.8% 


(6,1) 


15.8±0.1 


15.65 


0.93% 


(7,2) 


19.9±0.1 


20.11 


1.1% 


(7,3) 


20.8 ±0.1 


21.23 


2.1% 


(7,4) 


21.9±0.1 


22.42 


2.4% 


(8,1) 


22.5±0.1 


22.53 


0.14% 


(9,1) 


25.9±0.1 


25.97 


0.27% 


(9,2) 


26.6±0.1 


26.73 


0.5% 


(9,3) 


27.2±0.1 


27.55 


1.3% 


(9,5) 


30±0.1 


30.77 


2.6% 


(10,1) 


29.4±0.1 


29.41 


0.034% 


(10,2) 


29.9 ±0.1 


30.05 


0.52% 


(10,6) 


35±0.1 


35.93 


2.7% 


(20,1) 


73.78 ±0.01 73.83 


0.061% 


(20,4) 


73.79 ±0.01 74.01 


0.3% 


(20, 10) 


81.7±0.1 


82.57 


1.1% 


(20,20) 


109.7 ±0.1 


113.8 


3.7% 


(35,4) 


154±0.1 


154.1 


0.078% 


(35,10) 


158.2±0.1 


158.6 


0.27% 


(35,20) 


178.6±0.1 


180.4 


1% 



Superfluid State 

{N a ,N b ) E FND mc E A slda (error) 



(1,1) 


2.002 ±0 


2.302 


15% 


(2,2) 


5.051 ±0.009 


5.405 


7% 


(3,3) 


8.639 ±0.03 


8.939 


3.5% 


(4,4) 


12.573 ±0.03 


12.63 


0.48% 


(5,5) 


16.806 ±0.04 


16.19 


3.7% 


(6,6) 


21.278 ±0.05 


21.13 


0.69% 


(7,7) 


25 .923 ±0.05 


25.31 


2.4% 


(8,8) 


30.876 ±0.06 


30.49 


1.2% 


(9,9) 


35.971 ±0.07 


34.87 


3.1% 


(10,10) 


41.302±0.08 


40.54 


1.8% 


(11,11) 


46.889 ±0.09 


45 


4% 


(12,12) 


52.624 ±0.2 


51.23 


2.7% 


(13,13) 


58.545 ±0.1 8 


56.25 


3.9% 


(14,14) 


64.388 ±0.31 


62.52 


2.9% 


(15.15) 


70.927 ±0.3 


68.72 


3.1% 


(1,0) 


1.5 ±0.0 


1.5 


0% 


(2,1) 


4.281 ±0.004 


4.417 


3.2% 


(3,2) 


7.61 ±0.01 


7.602 


0.1% 


(4,3) 


11. 362 ±0.02 


11.31 


0.49% 


(7,6) 


24.787 ±0.09 


24.04 


3% 


(11,10) 


45.474 ±0.15 


43.98 


3.3% 


(15,14) 


69. 126 ±0.31 


62.55 


9.5% 



Table 9.2 Comparison between the ASLDA density functional as described in this section and the 
FN-DMC calculations 1 136 137 1 for a harmonically trapped unitary gas at zero temperature. The 
normal state energies are obtained by fixing A = in the functional: In the FN-DMC calculations, 
this is obtained by choosing a nodal ansatz without any pairing. In the case of small asymmetry, 
the resulting "normal states" may be a somewhat artificial construct as there is no clear way of 
preparing a physical system in this "normal state" when the ground state is superfluid. 



9.3.3.1 Trapped Systems 

The functional form of both the SLDA and ASLDA have been completely fixed by 
considering only homogeneous matter. Hence, a non-trivial test of the theory is to 
compare the energy of trapped systems with Monte Carlo calculations. This was first 



done for the SLDA in [ 127] and the results are shown in Fig. 9.12 Even for systems 
with only a few particles — which have large gradients — the agreement is within 10%. 
This rapidly improves to the percent level as one move to larger systems. 

The agreement is somewhat remarkable. In particular, we have included no gradi- 
ent corrections in the theory beyond the Kohn-Sham kinetic energy. These gradient 
corrections will contribute at some level, but in the present system the coefficients are 
extremely tiny (the leading gradient correction ~ (V«) 2 /« should give corrections 
that scale as E oc N 2 ^ for which there is no evidence in the Monte Carlo data). In 
any case, the agreement provides strong evidence that the SLDA captures the relevant 
energetics to provide a quantitative model of the unitary Fermi gas. 



45 



9.3 - Density Functional Theory for the Unitary Fermi Gas 



We should point out that the gradient terms in the SLDA are completely char- 
acterized by the kinetic terms. Thus, finite size effects are highly sensitive to the 
inverse effective mass parameter a. As mentioned in Sec. 9.3.2.1 the energy and 
gap can be fit with a = 1, but the resulting parametrization demonstrates a marked 
systematic deviation from the trap energies shown in Fig. 9.12 It is reassuring that 
the agreement is restored when the effective mass is chosen ( 9.95a| > to reproduce the 
quasiparticle spectrum. 



We have validated the ASLDA in a similar manner for trapped systems in Table 9.2 



Again, the agreement is at the few percent level in virtually all cases. In general, 
the formulation of the unitary DFT has a remarkable ability to capture the finite 
size effects in systems down to even a few particles [138], lending credence to the 
approximation of neglecting further gradients beyond the standard kinetic terms. 
This was somewhat anticipated since the kinetic terms completely describe finite 
size (shell) effects in the non-interacting system, but is non-trivial in the strongly 
interacting case of the unitary gas. 

Note that the BdG and SLDA functionals have also been considered in larger 
trapped systems M139II140I . 



9.3.3.2 Vortex Structure 

The first use of the SLDA was to determine the structure of superfluid vortices 11601 . In 
this work, two forms of SLDA (slightly different parameter values) were considered, 
and the solution for an axial symmetric vortex with unit circulation was found. 
The method of solution uses a technique that properly treats the infinite boundary 
conditions without truncating the physical space and introducing finite-size artifacts 
(see 111411 11421 for details). The profile for this vortex is shown in Fig. 9.13 In 



particular, it was predicted that the vortices should have a significant density depletion 
in the core — something that is not observed in the weak-coupling limit where pairing 
is exponentially suppressed. This predicted core depletion allows for the direct 
imaging of vortices in rotating trapped gasses 014311 . providing direct evidence for 
superfluidity in these systems. 



9.3.3.3 FFLO/LOFF 

The first application of the ASLDA was to consider the energetic stability of a 
Larkin-Ovchinnikov-Fulde-Ferrell (LOFF) 1146111471 polarized superfluid state 11251 . 
The density functional as constructed naturally supports a strong first-order phase 
transition between the fully paired superfluid state and the interacting normal state 



(dashed line in Fig. 9.14 1. We seed the functional with a periodic solution of the 
form shown in Fig. 9.15 with a node in the gap A (z). Allowing the system to relax 
to the optimal period L we find that this Larkin-Ovchinnikov type of solution has 
significantly lower energy than the competing pure and mixed phases over a large 
range of the phase diagram. 
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2 4,6 8 2 4 6 8 

k F P k p p 



Fig. 9.13 Density profile (left) and gap parameter (right) from 1 60 1 for a superfluid vortex in 
the symmetric n a = «/, unitary Fermi gas with unit circulation. The solid curve corresponds to a 
parametrization of the SLDA with no self-energy /3 = but including an effective mass correction. 
The dotted curve corresponds to a version with unit effective mass a = 1. The other two parameters 
were fixed to reproduce the best approximation to energies of the normal and superfluid states known 
at the time: Z,^ = 0-54 and ^sf — 0.44. The current parameter set l |9.95[ l should be preferred, but 
gives similar results. Note: The solid curve does not have the required currents to restore Galilean 
invariance (see Sec. |9.4.2fr , but the effect should be small here. Since the dotted curve has no 
effective mass correction, Galilean corrections are not required. 



This is a qualitatively new prediction of the ASLDA: such states are only meta- 
stable in the BdG. The effect of the self-energy corrections is to reduce the energy of 
these states to favor them over the homogeneous phases. It is interesting to note that 
the density contrast of these states is comparable to the density contrast in vortices 
(see Fig. 9.13 i. Such states have yet to be observed in experiments: this may be 
because the physical region in which the LO state is favoured exists only in a thin 
shell. Also, the one-dimensional structure discussed here will be unstable at any 
finite temperature [ 148 1 (see also 1 149 1) but might be stabilized in traps. The ground 
state will most likely be some sort of three-dimensional lattice structure (see for 
example [150]) which will likely require a fairly large physical volume to exists 
without significant frustration. The ideal situation would be a very flat trap tuned so 
that the LO region occupies a large physical space at the center of this trap, however, 
the construction of such traps presently poses some experimental difficulties that we 
hope will be overcome in the near future. 
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0.0 0.2 0.4 0.6 0.8 1.0 

x = n b /n a 

Fig. 9.14 The dimensionless convex function g(x) | 101 1 that defines the energy density S(n a ,ni,) = 
| (6tt 2 ) 2 / 3 [f! a ,g(x)] 5 ^ 3 as a function of the asymmetry x = «f,/n n (this plot is very similar to 
Fig. '|9.9) . The points with error-bars (blue online) are the Monte Carlo data from I132lll44ll64l . The 
fully-paired solution g(l) = (2^) 3//5 is indicated to the bottom right, and the recent MIT data (1331 
is shown (light x) for comparison. The phase separation discussed in 1 132, 144 64 101 1 is shown 
by the Maxwell construction (thin black dashed line) of the first-order transition. The LO state (thick 
red curve) has lower energy than all pure states and phase separations previously discussed. The 
Maxwell construction of the weakly first-order transition between the superfluid and LOFF phase is 
shown by the thick dashed line (red). 



9.4 Time-Dependent Superfluid Local Density Approximation 

9.4.1 Time-Dependent Equations for the Quasiparticle Wave 
Functions 



The equations for the time-dependent quasiparticle wave functions m„ ](T (r, t ) , v„ j(J (r , t ) 
have the time-dependent Bogoliubov-de Gennes form 



dt 





/ 


U b 




Vc, 




w 


V 



,+U a A 

h b + U b -A 
-A* -h*-U a 
A* -h%-U b J 





( U a \ 




Ub 




Va 


) 


w 



(9.100) 



For the sake of simplicity, we have dropped the arguments (r,f) for all functions 
in these equations. Note also that the external potentials U a (r,t) are real. The only 
difference with the static SLDA in the structure of h a (r,t) are the contributions arising 
from the variation of the current density correction to the kinetic energy density 
f (r,f ) ( 9.104|9.1 10 1, which are required by Galilean invariance to be discussed in 
Sec. 9.4.2 The chemical potentials jU fl ib , which can always be thought of as external 



constraints, are implicitly included in t/ CT (r,f). The chemical potentials can also be 
removed by a simple gauge transformation of the quasi-particle wave functions. It is 
straightforward to show that these equations conserve the total particle number for 
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-5 5 

z/k 



Fig. 9.15 A single LO period showing the spatial dependence of the pairing field A(z) (top) and the 
number densities of the majority (dotted) and minority (solid) species (bottom) at the values. Units 
are fixed so that = fl a — jllj is held fixed as it is for trapped systems. We normalize everything in 
terms of the density no = n a = interparticle spacing Iq = no -1 / 3 , and superfluid pairing gap Aq 
of the fully paired superfluid at the superfluid/LO transition point close to the center of the cloud. At 
the superfluid/LO transition, the character of the solution is that of widely spaced domain walls (see 
for example 11451 ). As one proceeds outward in the trap, the period and amplitude of the solution 
decreases until it is almost sinusoidal at the transition point to the interacting normal phase. 



arbitrary time-dependent external fields and also for arbitrary time variations of the 
coupling constant g(t). As expected however, in the presence of an external pairing 
field, particle number is not conserved: particles can be exchanged with the coupled 
system implied by the source of the external pairing field. 



9.4.2 Galilean Invariance 



The functionals as expressed in Sec. 9.3 are not manifestly covariant under Galilean 
transformation (a subset of the general coordinate invariance discussed in 01511 
which restrict the form of higher-order gradient terms). To restore this covariance, the 
currents currents j u (r) and jf,(r) described in (9.72 1 must be included. These vanish 



in time-reversal invariant ground states, but are crucial for discussing states that 
break time reversal and for the general time-dependent analysis. In nuclear physics 
Galilean covariance have been considered for quite some time II 1521 1 1531 1 1541 [1551 . 
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and the contribution of these currents is often essential for describing the properties 
of excited states. 

We start by expressing the Galilean invariance of the Lagrangian density for a 
single Fermi species (see 115111 for a more general discussion) 



> 2 \ 

\jf. (9.101) 



This is invariant under the following Galilean transformation: 

y/(x,f) -!>exp[-i(±m|v| 2 f + m\-x) /h] y/(x + vf,f)- (9.102) 

From this, we see that the currents and kinetic densities transform as 

j = I^Vl/A+h.c. -»j+mvn, (9.103a) 

t= 5 LVy/ r Vi//'->-T+v-j + im|v| 2 n. (9.103b) 

It follows directly that for a two-component system, the following combinations are 
Galilean invariant: 

|jcr| 2 jb ja (9.104) 



2m a n a m h n h m a n a 

We would like to separate out the center of mass motion from the intrinsic functional, 
so we introduce the total mass current and density: 

j+ = ja +jb, P+ = m a n a + m b n b . (9. 105) 

We may then write the functional in the following way: 

I . i"> 

Ij-f 



2p, 



■S. (9.106) 



The first term captures the energy of the center of mass motion and & describes the 
remaining intrinsic energy of the system, and should be strictly Galilean invariant. 

Excited states may be described by an extension of the DFT method commonly 
referred to as Time-Dependent Density Functional Theory (TDDFT) 1115611157111581 
159|. This theory describes the evolution of the one-body number density in the 
presence of an arbitrary one-body external field. As in the case of static DFT, one 
can prove an existence theorem I156II1571I1581 . This states that a functional exists 
from which one can determine the exact time-dependent number density for a given 
quantum system, and can be expressed in the form 
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h 2 



dtd 3 ri — £T ff (r,f)+£t/ ff (r,f)7i<r(r,f) 

[ ^ m a a 

+ ^[M a (r,0,f a (r,0,«6(r,0,T 6 (r,0,|v(r,?)| 2 ! g(r,?)]|. (9.107) 

Here a = a,b labels the two fermion species. The existence proof for superfluid 
systems is analogous to the proof for normal systems [ 1561 1157111581 . Here U a (r,t) 
are arbitrary time-dependent one-body external fields, which couple to the conserved 
number densities of the fermion species n a (r,t). These external fields can represent 
couplings to the laboratory environment, such as a trapping potential, which can be 
used to manipulate and study these systems. 

The direct coupling of an external gauge field to the electric charge and magnetic 
moments of the particles can also be incorporated in a straightforward manner, by 
the usual process of converting the global particle number symmetry to a local 
symmetry by invoking the principle of gauge invariance. We can also couple an 
arbitrary time-dependent external pairing field as well to represent interactions with 
another superfluid system brought into the proximity of the system under study. As 
mentioned above, this will violate the conservation of particle number as particles 
are now able to be exchanged with the other system. Finally, the last argument g(r, t) 
of the interaction term S represents the possibility of varying the coupling constants 
in space and time. In particular, as discussed in Sec. |9.2.1| by means of a Feshbach 
resonance an external magnetic field can be used to directly control the scattering 
length, providing yet another handle to manipulate and study these systems. 

In the functional 5? we have separated the kinetic energy /j 2 £ a T a (r,f)/2m from 
the interaction energy in order to disentangle the dependence on the reference frame. 
The interaction energy encoded in the functional § should be independent of the mo- 
tion of the system as a whole. By default, the properties of the ground states of a phys- 
ical system are typically discussed in the center of mass reference frame. When the 
system is excited by various external probes, inevitably currents appear. In the LDA it 
is natural to assume that the energy density separates into the kinetic energy of center 
of mass (which depends only on its local center of mass velocity and its corresponding 
mass) and the internal energy (which should not depend on the local center of mass 
velocity). The energy density S [n a {r,t),t a (r,t),n b (r,t),t b {r,t),\v(r,t)\ 2 ,g(r,t)] is 
the same as in the static SLDA, with the only difference that the dependence on the 



modified kinetic energy density f a (r,t) now includes the current densities (9.72b i to 
satisfy Galilean invariance 

ih 

Ja(r,0 = T E[Vv^(r,0<a('.0-v»,a(r,0Vv; ff (r ) 0] • (9.108) 

* n 

Here we will describe a slightly different philosophy in implementing the Galilean 
invariance than discussed at the beginning of this section, which leads to a somewhat 
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different definition of the modified kinetic energy densities f a (r,t) than those intro- 
duced in ( 9.104) > above. This ambiguity illustrates the freedom one has in introducing 



currents and using no other restriction except Galilean invariance. 

Upon boosting the system to a frame with a velocity V, the current density changes 
jcr(r,f ) — > j<j(r,f ) +mn(r,t)V. We introduce the velocity of the local center of mass 
frame (m a = my = m) 

V(r,.)-£§|, (9-109) 



where we have introduced the total current j + and density p + from (9.105 I. Conse- 
quently, the following combination of the kinetic energy density, current density and 
number density 

f a (r,t) = T ff (r,r) - j p (r,r) ■ V(r,r) + ^^( r »0V 2 (r,f) (9 nQ) 

renders the energy density locally Galilean invariant [127]. f CT (r,f ) is therefore the 
internal kinetic energy density in the local center of mass frame, which is different 



from the form of modified kinetic energy introduced in (9.1041. The difference 
between the two approaches to enforcing t he Gal ilean invariance amounts to terms 
proportional to \j b /m b n b -j a /m a n a \ 2 , see (9.104 1. 



It is worth noticing that because the Galilean invariance is built in, one of the 
famous relations in the Landau's Fermi liquid theory linking the effective mass of 
the quasiparticles with the p-wave interaction (denoted F\) is automatically satisfied 
(see JT2)). 

Note also that if terms arise of the form j fl (r,f) • j/^r,?), a new physical effect 
appears whereby the local velocity of one species depends also on the velocity 
of the other species. In other words, the inverse mass becomes a tensor in the spin 
("isospin") space. By including terms of the form |(j„(r.f) ■ Vn a (r,t))\ 2 , the effective 
mass becomes a tensor in real space. This was discussed in [160| in connection 
with the construction of the optimal local Schrodinger equation to represent a non- 
local equation. In particular, it seems that, in order to describe some rather subtle 
level orderings of the single-particle spectrum found in the a non-local Schrodinger 
equation, one needs a tensor effective mass in the local equations. This is also related 
to the discussion of superfluid mixtures, where it was observed long ago that one 
superfluid can drag the other one without any dissipation: the Andreev-Bashkin 
effect 11611 . Similar effects arise when one considers the terms induced by Galilean 



invariance (9.110* or j a (r,f) -j/,(r,f), when the presence of a current of one species 



induces a current of the other species. 



9.4.3 The Excitation of the Pairing Higgs Mode 

We shall illustrate the power of the Time-Dependence SLDA(TD-SLDA) by examin- 
ing the response of a superfluid unitary gas to the time variation of the scattering 
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length 1 162 1 . This problem has been studied extensively in the weak coupling regime 
whenk F \a\ < 1 and a < 0, see ifTBH [T64l [1651 fl66l fl67l [TBIfl [1691 [TTOl fT7T1fT72l 
[1731 El Q75] [1761 E21 UM EU [J80l HUl . The initial state of the system will be 
the ground state, and at subsequent times, the evolution will be adiabatic in the sense 
that no entropy production is allowed. To some extent this is a rather strong limitation 
of this time-dependent description of the quantum evolution, a restriction which can 
be lifted if one would consider a further extension of the formalism, the Stochastic 
TD-SLDA 1 182 1 which will not be discussed here. 




Fig. 9.16 The profile of the energy of a unitary Fermi gas as a function of the pairing gap with 
respect to the energy of the ground state. One would nai'ely expect that this system if released from 
a point almost at the tip of the "Mexican hat" will roll down along the radial direction, past the 
equilibrium value Aq = 0.5£f and oscillate indefinitely back and forth. 



Consider the following scenario 11621 : start with a homogeneous unitary Fermi 
gas in its ground state. At first slowly reduce the coupling constant y from its 
unitary value to a very small but still negative value. If this change is slow enough, 
then the system tracks the ground state into the ground state of the system with 
an exponentially small pairing gap. Now rapidly ramp the coupling y back to its 
unitary value and let the system evolve. This essentially looks at the evolution of 
an almost normal system with the unitary DFT. The behavior shown in Fig. 9.17 is 
rather surprising. 

Many approaches have been developed to describe the dynamics of a fermionic 
superfluid at or near T = including superfluid hydrodynamics, a Landau-Ginzburg 
or Gross-Pitaevskii (lg/gp) like description, and effective field theory, see 01831 
[IM[333[3l][l53][Il^ The common thread in all these approaches is 

the desire to identify a significantly smaller set of relevant degrees of freedom, and 
achieve an accurate description of the phenomena within a reduced framework. As 
a rule, when reducing the number of the degrees of freedom, one assumes that the 
system evolves through states where local equilibrium is maintained. In this instance, 
one would naively expect that the system dynamics are governed by an effective 
"Mexican hat" potential, Fig. 9.16 representing the energy of the system as a function 
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Fig. 9.17 The panels a,b and c display response of the homogeneous system to an initial 
switching time interval fo£f = 160, 10 and 160 and values of the gap corresponding to y s are 
J s /Y = 0.005 , 0.005 and 0.5 respectively, where y is the coupling constant controlling the magnitude 
of the pairing gap and Aq rs 0.5£f is the gap equilibrium value, both at unitarity. The panels a and 
b show that when the system is released from the neighborhood of the tip of the "Mexican hat" 
potential the pairing gap oscillates back and forth, but never past the equilibrium value Aq. At the 
same time the system will rotate around the origin as the phase of the pairing field (not shown here) 
will monotonically evolve in time as well. However, when the system is released from an initial 
position closer to the minimum at Aq the oscillation is damped, A (/) = A^ +A s'm(2A„t + <j>)/ \/A x t, 
with a mean frequency 2A X and around a value smaller than the equilibrium A x < Aq, a behavior 
which was first studied in 1 163| in the BCS limit, when the coupling is weak. 

of the complex pairing field. The system is brought adiabatically from the minimum 
of the potential to almost the "tip of the Mexican hat", and released with zero initial 
"velocity". The naive picture is that the system will "roll" down along the radial 
direction accelerating until it reaches the minimum of the potential, pass through the 
minimum, and oscillate back and forth along the "radial" direction without damping. 

One might also inspect the LG/GP description of the dynamics of the system using 
the nonlinear Schrodinger equation 

m ^ = _^o + ^ (lnri0 , 2)nr , . (9 . m) 

Since there are no spatial gradients in this system (we have changed the coupling in 
a uniform manner so as not to break the translational invariance), only the second 
term on the right hand side of this equation contributes, and the solution is a simple 
monotonic evolution of the condensate phase f (r,f) in time and the magnitude of 
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"wave function" f^r,?) remains constant. In lieu of a better simple alternative, many 
authors have used this approach to characterize dynamics of Fermi superfluids at 
essentially zero temperatures, even though the LG/GP description is only justified 
near the critical temperature. 

Another approach is to use the zero-temperature limit of Landau's two fluid 
hydrodynamics, which reduces to the following two equations at zero temperature 

{m\ 2 (r t) 1 
2 { m ' +fi[w(r,0]}=0. 

(9.112) 

Here v(r,f) is the hydrodynamic velocity and ji[n(r,t)] is the local thermodynamic 
potential. Since there are no spatial gradients, these two equation simply predict that 
the number density will remain constant and nothing else will happen. 




Fig. 9.18 Panels a and b display the instantaneous occupation probabilities of the mode shown 
in upper panel of Fig. |9. 17| corresponding to times t > when the pairing field is at its minimum 
and maximum values respectively with circles joined by a solid (blue with circles) line. With (red) 
dots we plotted the equilibrium occupation probabilities corresponding to the same instantaneous 
values of the pairing gap. In panels c and d we show the maximum and minimum values of the 
oscillating pairing field and the corresponding excitation energy as a function of the frequency of 
the Higgs-like modes, see Fig. |9.17| q and b. 



Apart from the fact that the number density will remain constant and spatially 
uniform, these three different naive pictures lead to drastically different predictions. 
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The actual time evolution of the system, shown in Fig. |9.17| is qualitatively different 
from any individual picture, but demonstrates a combination of the expected features. 
The pairing gap does increase from almost zero towards the equilibrium value, and it 
oscillates, but it never crosses the minimum equilibrium value Aq of the "Mexican 
hat" potential. At the same time, the phase of the pairing gap increases monotonically 
in time and the number density is constant. 

By preparing the initial state slightly differently one can excite different types of 
these modes that have been dubbed "Higgs" modes of the pairing field. One can vary 
the upper and the lower values between which the pairing field oscillates, and also 
adjust the period of these oscillations. It is remarkable, however, that the frequencies 
of these modes are always smaller than 2Aq, where Aq is the equilibrium value of the 
pairing gap at unitarity, even though the excitation energy is large. These are indeed 
very collective excitations of the system, of extremely low frequency, but with an 
excitation energy per particle significantly less than pairing gap. 

It is still an unresolved question of how these modes will eventually decay and 
how the system will thermalize. It is also instructive to examine the time dependent 
population of the various single-particle momentum states of these collective modes 
as shown in Fig. 9.18 When the value of the pairing gap is very small the occupation 
probabilities are essentially those of a system in equilibrium. However, when the 
system reaches a pairing gap essentially equal to the equilibrium value Aq, the 
occupation probabilities are clearly very different from those in the ground state, 
which clearly points to a non-equilibrium state. This aspect should clarify why neither 
LG/GP nor quantum hydrodynamics are valid as both assume local equilibrium is 
maintained. 



9.4.4 Generation and Dynamics of Vortices 

A number of results concerning the generation and dynamics of vortices in a unitary 
Fermi gas by an external time-dependent perturbation can be found at [ 1 88 1 . As far as 
we are aware, this problem has been studied in one paper for a pure 2D systems 1 1 89 1 . 
As in the previous example, we do not yet consider entropy production in these 
simulations. 

In order to illustrate further the power of the TD-SLDA as well as the limitations 
of traditional approaches such as superfluid hydrodynamics or a LG/GP analysis, 
we now consider the quantum dynamics of a stirred unitary Fermi gas [ 1 88 1 . We 
start with the gas in its ground state in a cylindrical trap, uniform and with periodic 
boundary conditions in the third spatial direction. We then subject the system to a 
time-dependent external stirring field which breaks the cylindrical symmetry. When 
implemented numerically H 1901 , if one places the system on a spatial lattice with 
N s spatial lattice points in one direction, one can show that the size of the problem 
scales as oc jy~. When the limitation of spatial homogeneity in the z-direction is lifted 
the size of the problem scales as oc Nf, which as a rule requires an implementation 
on the largest leadership class supercomputers available. For example, if N s = 50 an 
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efficient solution of the TD-SLDA equations becomes possible only on the JaguarPF 
Cray XT50which we are currently utilizing to its full extent. 

When homogeneity along the z-direction is enforced, the quasiparticle wave 
functions have the structure (u n (x,y,t) exp(ikz) ,v n (x,y,t) exp(ikz)) while self -energy 
U(x,y,t) and the pairing potential A(x,y,t) are translationally invariant along z. We 
adiabatically introduce a vertical rod into this "soup can" and stir the gas with a 
constant angular velocity. One can vary both the stirring radius R and stirring angular 
frequency co to control the speed v ro d = Rco of the rod. 

One might expect that if v ro( j -C v c , where v c is the critical velocity of a unitary 
Fermi gas, then the system will return to its initial state after the stirring is turned 
off. However, if v ro d > v c , then one might destroy the superfiuid order, resulting in a 
normal Fermi gas. If v\ < v ro d < v c , where vi is some minimal stirring velocity, one 
expects that vortices will be created. Unfortunately, none of the simple theories can 
shed much light on the outcome: superfiuid hydrodynamics cannot describe quantum 
vortices as there is no intrinsic quantization or Planck's constant in its formulation: 
vortex quantization must be imposed by hand, and there is nothing in principle to 
prevent decay of a quantized vortex into two fractionally quantized vortices. The 
time dependent LG/GP approach will also fail to describe the normal state and the 
transition from superfiuid to normal state, as it is formulated explicitly in terms of 
the order parameter alone, which vanishes in the normal state. Thus, it seems that 
the only viable solution is to forgo a reduction in the degrees of freedom and deal 
directly with the qu as ip articles included in the DFT. 

A unitary Fermi gas is a special system in quite a number of ways: in particular, 
it appears to have the highest critical velocity of all known superfluids ||191||192l . 
On the BCS side of the Feshbach resonance, if stirred fast enough, the system can 
loose superfluidity by the breaking of the Cooper pairs v qp — mia(Ei c ,k), On the BEC 
side of the Feshbach resonance, the dominant mechanism for the loss of superfluidity 
is the excitation of the Anderson-Bogoliubov sound modes c = /3. In the 

a unitary Fermi gas, these two different critical velocities appear to be essentially 
identical, and exactly at unitarity one obtains 




v c = rnin(c, v qp ) = y F min^y|,yay(j3-^) 2 + Tj2 + (j3-^)j sa 0.365v F . 
Since the amount of information one extracts in a TD-SLDA simulation of this type is 



very large, it is not sufficient to display only a few pictures such as those in Figs. 9.19 
and |9. 20 1 We invite the interested reader to explore some of the movies made of these 
simulations 11881 . We shall comment here only on a few selected aspects of these 
results, most of which will be prepared for a publication at a later time. 

Our expectation that, under gentle stirring, the unitary Fermi gas will return to 
its initial superfiuid state is supported by the simulations [ 188 1, and is in line with 
how one would expect a superfiuid to respond to such an external probe. The other 
expectation, that vigorous stirring can destroy the superfiuid order is also confirmed. 



' JaguarPF is a Cray XT5 supercomputer with 224,256 processing cores, see http://ww.nccs.gov 
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Fig. 9.19 The contour density profiles of a unitary Fermi gas in a cylindrical container, stirred with 
a uniformly rotating rod, which is inserted and extracted adiabatically from the system. The position 
(and intensity) of the rod can be inferred as the deepest density depletion in the system, and it is 
actually visible only in panels 2-4. Initialy the gas shows an almost uniform density distribution, 
and subsequently it is gathered almost entirely in front of the stirring rod. The magnitude of the 
density is in units of the unperturbed central initial density of the cloud and the colorbar on the right 
decodes the meaning of each color used. By the end of the simulation there are 13 vortices forming 
an almost perfect triangular Abrikosov lattice in this confined geometry. 

Within TD-SLDA, the dynamic generation of vortices as well as formation of the 
celebrated Abrikosov vortex lattice are also readily demonstrated. By varying shape, 
the number and the stirring velocity we generated a plethora of quantized vortices in 
this "soup can" of superfluid unitary Fermi gas H188II . 

While we expected to generate a relatively small number of vortices at low stirring 
velocity, and that the number of vortices will increase with more vigorous stirring, 
many of the features of dynamic vortex generation are quite surprising. The fact that 
this system is compressible results in surprisingly large time-dependent variations 
of the local number density. Often the entire mass of the system is gathered in 
front of the stirrer, leaving little matter behind it: The gas can occupy less than half 
of the available volume, even though the volume excluded by the stirrer is quite 
small. It also comes as a great surprise that the system does not loose quantum 
coherence under such a violent perturbation. Moreover, it organizes itself in an 
almost perfect vortex lattice after the stirring is turned off. Even more surprising is 
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Time [1/e F ]: 547 Time [1/e p ]: 730 Time [1/e p ]: 912 




Fig. 9.20 The corresponding contour profiles of the pairing field |4(x,v,/)| a unitary Fermi gas in a 
cylindrical container, stirred with a uniformly rotating rod. A plot (not shown here) of the phase 
of the pairing field arg A (x,y,t) shows that as one circles a vortex core the phase changes by In, 
thus each of these vortices carries exactly h/2 units of angular momentum per particle and both the 
number (normal) density and the pairing field are significantly depleted in the core of the vortex 
(60). 



that the system remains superfluid, even when stirred at supercritical speeds! We 
have observed that the system forms a vortex lattice even if stirred with speeds up to 
v s = 0.60v/r > V c ~ 0.365vf (see the case of 7 vortices with a large radius stirrer at 
1 188 1). We attribute this behavior to the fact that an increased density of the cloud 
during the stirring process corresponds to an increased local critical velocity, since 
the local Fermi velocity increases as well accordingly. 

These two cases of exciting and monitoring the unitary Fermi gas by two drasti- 
cally different methods illustrate both the power and flexibility of this framework, as 
well as the richness of the phenomena waiting to be fully explored. One potential 
topic to be explored by these techniques that has mesmerized the low temperature 
community during the past few decades is quantum turbulence 111931 11941 11951 . 
Hopefully this can also be replicated in experiments with cold atomic fermionic 
gases. Due to the complexity of the full 3D time-dependent Bogoliubov-de Gennes 
equations, this aspect has never been theoretically addressed for fermionic systems. 
The TD-SLDA appears as a framework of choice in this respect. In particular one 
can address on a fully microscopic basis for vortex reconnection dynamics, which 



59 



9.5 - Concluding Remarks 



is likely at the heart of quantum turbulence at zero temperature, where dissipative 
processes are greatly inhibited. 



9.5 Concluding Remarks 

We have reviewed here three methods to describe the properties of many-body 
systems starting from the bare Hamiltonian, and building a practical framework for 
studying nontrivial properties of mesoscopic systems and quantum dynamics. 

The QMC method is particularly suited to calculate properties of the homogeneous 
phase of matter in an unbiased fashion. It also can be used for inhomogeneous 
systems, but is limited by system size and can not handle large number of Fermions. 
It is also generally plagued by the infamous sign problem (except in exceptionally 
symmetric contexts) and so far has not been used to describe systems in the time 
dependent domain. 

The complimentary approach of density functional theory (DFT) through the use 
of the SLDA and ASLDA can be applied to extend these results to mesoscopic systems 
with larger number of particles and a wide variety of geometries. The time dependent 
TDDFT (TD-SLDA) extension brings these techniques to bear on time-dependent 
quantum dynamics. The main difficulty with the DFT is that there is no well defined 
procedure to construct the functional. However, in the particular case of a unitary 
Fermi gas, the form of the SLDA and ASLDA functionals is sufficiently restricted by 
dimensional analysis, QMC results, and Galilean invariance as to be able to make 
testable predictions with relatively small uncertainties. This has been validated with 
both ab initio theoretical and experimental results. 

The next step is to use such tools to make predictions about the properties of 
the unitary gas under various conditions: by changing the geometry and even the 
Hamiltonian as a function of time, by probing the system with a variety of external 
probes and exciting a plethora of modes — both linear and nonlinear — and by studying 
both the equilibrium and non-equilibrium dynamics. We have illustrated a few of 
these applications, but it is clear that we have barely scratched the surface of this 
subject. 

The Fermi gas in the unitary regime proves to be an extraordinarily rich physical 
system to study, not only because one can both theoretically and experimentally 
address many of its properties with both precision and accuracy, but because it has 
so many truly unexpected phases and dynamical phenomena. 

Many fascinating features of this systems are still waiting to be revealed in experi- 
ments in their full glory, including: the pseudogap phase, the supersolid LOFF phase, 
/7-wave superfluidity 1 1961 I197L the Higgs mode of the pairing field, the behavior and 
response to various spatial and time varying trapping fields and probes, the dynamics 
of vortices which opens a window to quantum turbulence both theoretically and 
experimentally, the existence of supercritical superflow, and its kinetic properties — in 
particular its viscosity. One can safely state that the most extraordinary features of 
the unitary gas are still waiting to be demonstrated. 
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Perhaps the most captivating part of the story of the unitary Fermi gas is that it 
provides a link to an abundance of widespread fields of physics, from optics and 
atomic physics, to condensed matter physics, to nuclear physics and the physics of 
neutron stars, color superconductivity in QCD and dark matter, relativistic heavy ion 
collisions, and the Ads/CFT approaches in quantum field theory. 
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9.6 Appendix 

9.6.1 Formal Description of the DFT 

Here we present a somewhat formal derivation of the variational property of the 
Kohn-Sham equations. Consider a general free-energy functional of the following 
form 

F=E(n Al n B ,---) + TTr(plnp) (9.113) 
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where 

« A =Tr(pA T ), 
« B =Tr(pB r ), 



are the various densities, anomalous densities, etc. expressed linearly in terms of 
the one-body density matrix p. By varying the functional with respect to the density 
matrix p subject to the appropriate constraints on density matrix form (discussed in 
section |9~.6.1.1| >, one obtains a solution of the form 

p=/ j8 (H[p]) (9.114) 

where fg (E) is the appropriate thermal distribution for the particles of interest, and 
H is a single-particle Hamiltonian that depends on p : 

, , dE dE 
H[p = — A + — B + • • • . (9.115) 

The typical Kohn-Sham equations follow by diagonalizing the self-consistency 
condition \9.\ 14| l with a set of normalized Kohn-Sham eigenfunctions of definite 
energy: 

U\n)=E n \n). (9.116) 

The density matrix is diagonal in this basis and expressed in terms of the appropriate 
distribution functions fp (E): 

p=Y,fp(E n )\n)(n\. (9.117) 

n 

All of the functionals considered in this chapter may be expressed in this form. Once 
the appropriate matrix structures A, B etc. are described, the form of the Kohn-Sham 
equations and potentials follows directly from these expressions. 



9.6.1.1 Fermions 

The only remaining complication is to impose the appropriate constraint on the 
density matrix p. This ensures that the appropriate statistics of the particles is 
enforces. As we shall be interested in Fermions, the relevant constraint on the density 
matrix (dictated by the canonical commutation relationships) is 

p + Cp r C = l (9.118) 

where C = C r is the charge conjugation matrix: 

C|vO = IV>*- (9-119) 
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This follows from the anti-commutation relationship for fermions and is discussed 



further in the Section 9.6.2 The constrained minimization of the functional F(p) 
results in the standard Fermi distributiorQ 



j3(H[p]-CH r [p]C) 



1 



(9.120) 



which is the fermionic form of the self-consistency condition (9.114i for the density 



matrix p. In practise, one does not iterate the entire density matrix. Instead, one 
stores only the densities ha, «b, etc. Through ( |9.115| l, these define the Kohn-Sham 
Hamiltonian H, which is then diagonalized to form the new density matrix and 
finally the new densities. If, for example, symmetries allow the Hamiltonian H to be 
block diagonalized, then one can construct and accumulate the densities in parallel 
over each block. Finally, the densities represent far fewer parameters than the full 
density matrix. Thus, more sophisticated root-finding techniques such as Broyden's 
method |199| may be efficiently employed: Applying these techniques to the full 
density matrix would be significantly more expensive. 



9.6.2 Single Particle Hamiltonian 



It is convenient to express these concepts in the language of second quantization. The 
Hamiltonian will appear as a quadratic operator of the form 



(9.121) 



where f has several components and s is a matrix. The factor of 1/2 accounts for 
the double counting to be discussed below. For a two component system, the most 
general f that allows for all possible pairings has four components: 



f = 



b 



(9.122) 



In terms of components of the wavefunction, we will write ,3V s y = Ey/ where: 



u b 

V a 

w 



(9.123) 



7 Formally, this constraint can be implemented using a Lagrange multiplier, but it is much easier to 
see the results by letting p = l/2 + x — Cx r C where x is unconstrained, and then performing the 
variation with respect to x. 
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The naming of these components is conventional (see for example [117]) and the 
functions u and v are typically called "Coherence Factors". Note that the conven- 
tion is that v* b (r,t) are the wavefunctions of the particles. In this formulation the 
Hamiltonian has the form presented in ( |9.100| i: 



(h a + U a A \ 

h b + U b -A 

-A* -h*-U a 
V A* -h* b -U b J 



(9.124) 



9.6.2.1 Four-component Formalism 

We shall start with this full four-component formalism but soon utilized a reduction: 
If the superfluid pairing A ~ (ab) channel is attractive, then often the "Fock" channel 
is repulsive so we can take (a'b) — 0. In combination with the double-counting 
discussed below, this will allow us to fully express the system in terms of two 
components. 

The four-component formalism double counts the degrees of freedom: f contains 
both a and a ' . This degeneracy is described in terms of the charge conjugation matrix 

Y.: 

f f = "iff where ^ = (l o) " (9.125) 

The operator f will satisfy the single-particle Shrodinger equations 

H S \W)=E\W) (9.126) 

where the Hamiltonian H s = f ^ Jtff ' s *¥ can be chosen to satisfy (the sign implements 
Fermi statistics) 

c gM' T s ^'=-M' s . (9.127) 

In this form, the charge conjugation symmetry ensures that the eigenstates will appear 
in ±E pairs [^Keeping only one set of pairs will ensure that we do not double count. 
Using this symmetry, we can formally diagonalize the Hamiltonian by a unitary 
transformation °k such that: 

&&& = \fc °J] (9.128) 



2 -e / 

where E = diag (£■,-) is diagonal. The columns of the matrix ^ are the (ortho)normalized 
wave-functions and describe the "coherence" factors. To determine the correct expres- 



8 Suppose .ffi's'ty = Elff. Applying 1 9. 127 , using ^ 2 = 1, and taking the transpose imply that 
\l/ Tc ^ T Jf j = — £i// r< if r . Since left and right eigenvalues are the same, this implies that there is 
some other state such that ,Jrff s \jr = — Bijf. For Hermitian Hamiltonians, Jff s = 3v[, hence, the other 
state can be directly constructed as if/ = c €\y* . 
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sions for the densities in terms of the wavefunctions we form them in the diagonal 
basis and then transform back to the original basis using a i/ . 

Despite this formal degeneracy of eigenstates, we are not aware of a general 
technique to block-diagonalize the original Hamiltonian in the presence of non-zero 
terms of the form (a'b), though perhaps the symmetry might be incorporated into 
the eigensolver. 



9.6.2.2 Two-component Formalism 

If {a'b) — 0, however, then the Hamiltonian is naturally block diagonal: 

= \ _rt), H ' = V tH * V + const . 

and one may consider only a single block in terms of the reduced set of operators 

(9.130) 



(9.129) 



a 



This directly avoids any double counting issues. This system may be diagonalized: 

H,U = UE. (9.131) 

The matrix U defines the single "quasi"-particle operators <j) as linear combination of 
the physical particle operators contained in y/: 



4> = uV- 

The Hamiltonian is diagonal in this basis 
and hence expectation values may be directly expressed 
((j>^) = e p (E) = 



(9.132) 



(9.133) 



V 



where 1 — 0^ (E) = fp (E) is the appropriate distribution function: For fermions we 
have 

W E)= uT*e- (9 ' 134) 

At T = this reduces to 6q(E) = 9(E) and is equivalent to the zero-temperature 
property that negative energy states are filled while positive energy states are empty. 
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This may be simply transformed back into the original densities (on the diagonal) 
and anomalous densities (off-diagonal): 

Fermi statistics demands F +F + = 1 but we may have to relax this requirement 
somewhat in order to regulate the theory in terms of an energy cutoff d c (E). The 
columns of U„ of U correspond to the single-particle "wavefunctions" for the state 
with energy E n . We partition these into two components sometimes referred to as 
"coherence factors" 

u « =(;;)• (9 - i35) 

The unitarity of U imposes the conditions that 

u m u « + V m V n = 5 mn , (9.1 36a) 

£u„i4=£v„vj = l, (9.136b) 

n n 

£u n vj = £v„i4 = 0. (9.136c) 

n n 

From this we may read off the expressions for the densities 

n a = (a f a)=Y J <u T „9 li (-E n ), (9.137a) 

n 

n h = (tfb) =Y,Vnvl9p(E n ), (9.137b) 

n 

v = {ab)=Y j yy n y%{E n ) (9.137c) 

n 

= -Y^» n y%(-E n ) (9.137d) 

n 

= Zu n vl 9p{En) - 9p{ - En) . (9.137e) 

n Z 

The last form for v must be used if the regulator is implemented such that 9 C (E) + 
d c (—E) 7^ 1, in particular, if 6 C (E) = for \E\ > E c . Note that these expressions are 
basis independent, e.g. in position space: 

» a (r,r') =J>„(r)*«„ (r'f dp (-£„). (9.138) 

n 

The energy E n here is the energy determined by solving these equations and will 
contain both positive and negative energies. 
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